#include "Func.hpp"
#include "KdTree.hpp"
#include "cmsh.h"

#include <iostream>
// #include <cmsh.h>
// #include <mpi.h>

using namespace std;

/*
                                        6
                z ▲ +------------------------------------>+
                |  /▲                                   /\▲
                7 / |                                8 /  | 
                |/  | 2           5                   /   | 4
                +---+------------------------------->+    |
                ▲   |                                ▲    |
                |   +--------------------------------+--->+
              1 |  /\              l2              3 |   /\
                | / l3                               |  / l4
                |/                                   | /
                +------------------------------------>+---------> y
            （0, 0）               l1

                    ‌II‌                 2                   IV
                    +--------------------------------+--->+
                   /\                                    /\
                3 /            Surface{1}               / 4  
                 /                                     /
                +------------------------------------>+
                I                   1                  ‌III‌  
*/
template <typename datatype>
int GridRemesh(const CornerGrid<false, datatype, GridCell<GridPoint<datatype>, datatype, 8>>& basegrid,
                FinElemGrid<GridCell<int, datatype, 4>, datatype, 4>& dstgrid, 
                double lc, double xExtensionFactor, double yExtensionFactor, double zExtensionUpper, double zExtensionLower)
{
    int num_surf = static_cast<int>(basegrid.GetNx() * 1.0);
    if (num_surf < 1) {
        num_surf = 1;
    }
    
    int point_idx = 0, line_idx = 0, surface_idx = 0, volume_idx = 0;
    int LastPoint_1 = 0, LastPoint_2 = 0, LastPoint_3 = 0, LastPoint_4 = 0;
    int LastPoint_5 = 0, LastPoint_6 = 0, LastPoint_7 = 0, LastPoint_8 = 0, LastPoint_9 = 0, LastPoint_10 = 0, 
        LastPoint_11 = 0, LastPoint_12 = 0, LastPoint_13 = 0, LastPoint_14 = 0, LastPoint_15 = 0, LastPoint_16 = 0;
    int LastLine_1 = 0, LastLine_2 = 0, LastLine_3 = 0, LastLine_4 = 0,
        LastLine_x1 = 0, LastLine_x2 = 0, LastLine_x3 = 0, LastLine_x4 = 0;
    int LastLeftLine_1 = 0, LastLeftLine_2 = 0, LastRightLine_1 = 0, LastRightLine_2 = 0;
    int LastFrontSurf = 0, LastBackSurf = 0;    
    int NX = basegrid.GetNx(), NY = basegrid.GetNy(), NZ = basegrid.GetNz();

    std::vector<int> EnterLastPoint, LastSurface;
    std::vector<int> LastLeftSurf, LastRightSurf;
    std::vector<int> compoundLines_y1, compoundLines_y3, compoundLines_y2, compoundLines_y4;
    std::vector<int> LeftXpoint[2], RightXpoint[2], LeftXpoint_ext[2], RightXpoint_ext[2];
    std::vector<int> LastxLine_1, LastxLine_2, LastxLine_1_ext, LastxLine_2_ext;
    std::vector<int> compoundLines_z1, compoundLines_z2, compoundLines_z3, compoundLines_z4;

    std::vector<int> cellzone_id[6];
    std::vector<int> Interface_surf[NZ + 3];
    std::vector<int> Interface_line[NZ + 3];

    datatype X_max, X_min, Y_max, Y_min, Z_max, Z_min;
    datatype X_ext_max, X_ext_min, Y_ext_max, Y_ext_min, Z_ext_max, Z_ext_min;

    X_max = basegrid.GetXmax(), X_min = basegrid.GetXmin();
    Y_max = basegrid.GetYmax(), Y_min = basegrid.GetYmin();
    Z_max = basegrid.GetZmax(), Z_min = basegrid.GetZmin();

    X_ext_max = X_max + xExtensionFactor * (X_max - X_min);
    X_ext_min = X_min - xExtensionFactor * (X_max - X_min);
    Y_ext_max = Y_max + yExtensionFactor * (Y_max - Y_min);
    Y_ext_min = Y_min - yExtensionFactor * (Y_max - Y_min);
    Z_ext_max = -zExtensionLower;
    Z_ext_min = -zExtensionUpper;

    cmsh::initialize();

    auto CellArray = basegrid.GetCellArray();
    auto FaultArray = basegrid.GetFaultArray();
    
    cmsh::clear();
    cmsh::model::add(basegrid.GetGridName());

    cmsh::option::setNumber("Mesh.Algorithm", 5);
    // cmsh::option::setNumber("Mesh.Optimize", 1);
    // cmsh::option::setNumber("Mesh.OptimizeNetgen", 1); // 使用 Netgen 进一步优化
    cmsh::option::setNumber("Mesh.CharacteristicLengthMin", lc);
    cmsh::option::setNumber("Mesh.CharacteristicLengthMax", 1.333 * lc);

    /* 1、Divide geological interfaces */
    // 下层
    for (int i = 0; i < 2; i++) {
        LeftXpoint[i].resize(num_surf + 1);
        RightXpoint[i].resize(num_surf + 1);
        LeftXpoint_ext[i].resize(num_surf + 3);
        RightXpoint_ext[i].resize(num_surf + 3);
    }
    LastxLine_1.resize(num_surf);
    LastxLine_2.resize(num_surf);
    LastxLine_1_ext.resize(num_surf + 2);
    LastxLine_2_ext.resize(num_surf + 2);

    int bcIdx_1 = 0 * NX         + 0;
    int bcIdx_2 = 0 * NX         + NX - 1;
    int bcIdx_3 = (NY - 1) * NX + 0;
    int bcIdx_4 = (NY - 1) * NX + NX - 1;
    cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[0]._x, CellArray[bcIdx_1]._points[0]._y, Z_ext_min, lc, 1);
    cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[1]._x, CellArray[bcIdx_2]._points[1]._y, Z_ext_min, lc, 2);
    cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[2]._x, CellArray[bcIdx_3]._points[2]._y, Z_ext_min, lc, 3);
    cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[3]._x, CellArray[bcIdx_4]._points[3]._y, Z_ext_min, lc, 4);

    cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[0]._x, Y_ext_max, Z_ext_min, lc, 5);
    cmsh::model::geo::addPoint(X_ext_min, CellArray[bcIdx_1]._points[0]._y, Z_ext_min, lc, 6);
    cmsh::model::geo::addPoint(X_ext_min, Y_ext_max, Z_ext_min, lc, 7);

    cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[1]._x, Y_ext_max, Z_ext_min, lc, 8);
    cmsh::model::geo::addPoint(X_ext_max, CellArray[bcIdx_2]._points[1]._y, Z_ext_min, lc, 9);
    cmsh::model::geo::addPoint(X_ext_max, Y_ext_max, Z_ext_min, lc, 10);

    cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[2]._x, Y_ext_min, Z_ext_min, lc, 11);
    cmsh::model::geo::addPoint(X_ext_min, CellArray[bcIdx_3]._points[2]._y, Z_ext_min, lc, 12);
    cmsh::model::geo::addPoint(X_ext_min, Y_ext_min, Z_ext_min, lc, 13);

    cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[3]._x, Y_ext_min, Z_ext_min, lc, 14);
    cmsh::model::geo::addPoint(X_ext_max, CellArray[bcIdx_4]._points[3]._y, Z_ext_min, lc, 15);
    cmsh::model::geo::addPoint(X_ext_max, Y_ext_min, Z_ext_min, lc, 16);

    compoundLines_z1.push_back(1);
    compoundLines_z2.push_back(2);
    compoundLines_z3.push_back(3);
    compoundLines_z4.push_back(4);

    compoundLines_y1.push_back(1);
    compoundLines_y3.push_back(2);
    LeftXpoint[0][0] = 1;
    RightXpoint[0][0] = 3;

    LeftXpoint_ext[0][0] = 7;
    LeftXpoint_ext[0][1] = 5;

    RightXpoint_ext[0][0] = 13;
    RightXpoint_ext[0][1] = 11;

    int LastxPoint_1 = 1, LastxPoint_2 = 3;
    int LastxPoint_1_ext = 5, LastxPoint_2_ext = 11;

    point_idx += 16;

    for (int j = 0; j < NY - 1; j++) {
        // i = 0
        bcIdx_1 = j * NX + 0; 
        // i = NX - 1
        bcIdx_2 = j * NX + NX - 1;

        cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[2]._x, CellArray[bcIdx_1]._points[2]._y, Z_ext_min, lc, point_idx + 1);
        compoundLines_y1.push_back(point_idx + 1);
        cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[3]._x, CellArray[bcIdx_2]._points[3]._y, Z_ext_min, lc, point_idx + 2);
        compoundLines_y3.push_back(point_idx + 2);

        point_idx += 2;
    }

    compoundLines_y1.push_back(3);
    compoundLines_y3.push_back(4);

    cmsh::model::geo::addBSpline(compoundLines_y1, 1);
    cmsh::model::geo::addBSpline(compoundLines_y3, 2);

    cmsh::model::geo::addLine(6, 12, 3);
    cmsh::model::geo::addLine(9, 15, 4);

    LastLine_3 = 3, LastLine_4 = 4;

    line_idx += 4;

    cmsh::model::geo::addLine(5, 7, 5);
    LastxLine_1_ext[0] = 5;
    cmsh::model::geo::addLine(1, 5, 6);
    cmsh::model::geo::addLine(6, 7, 7);
    LastLeftLine_1 = 7;

    cmsh::model::geo::addLine(11, 13, 8);
    LastxLine_2_ext[0] = 8;
    cmsh::model::geo::addLine(3, 11, 9);
    cmsh::model::geo::addLine(12, 13, 10);
    LastRightLine_1 = 10;

    line_idx += 6;

    int LastyLine = 1, LastLefty = 6, LastRighty = 9;

    for (int i = 0; i < num_surf - 1; i ++) {
        // j = 0
        bcIdx_3 = 0 * NX + i;
        // j = NY - 1
        bcIdx_4 = (NY - 1) * NX + i;

        cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[1]._x, CellArray[bcIdx_3]._points[1]._y, Z_ext_min, lc, point_idx + 1);
        LeftXpoint[0][i + 1] = point_idx + 1;
        cmsh::model::geo::addLine(LastxPoint_1, point_idx + 1, line_idx + 1);
        LastxLine_1[i] = line_idx + 1;

        cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[3]._x, CellArray[bcIdx_4]._points[3]._y, Z_ext_min, lc, point_idx + 2);
        RightXpoint[0][i + 1] = point_idx + 2;
        cmsh::model::geo::addLine(LastxPoint_2, point_idx + 2, line_idx + 2);
        LastxLine_2[i] = line_idx + 2;

        cmsh::model::geo::addLine(point_idx + 1, point_idx + 2, line_idx + 3);

        // left_ext
        cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[1]._x, Y_ext_max, Z_ext_min, lc, point_idx + 3);
        LeftXpoint_ext[0][i + 2] = point_idx + 3;
        cmsh::model::geo::addLine(LastxPoint_1_ext, point_idx + 3, line_idx + 4);
        LastxLine_1_ext[i + 1] = line_idx + 4;

        cmsh::model::geo::addLine(point_idx + 1, point_idx + 3, line_idx + 5);

        // right_ext
        cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[3]._x, Y_ext_min, Z_ext_min, lc, point_idx + 4);
        RightXpoint_ext[0][i + 2] = point_idx + 4;
        cmsh::model::geo::addLine(LastxPoint_2_ext, point_idx + 4, line_idx + 6);
        LastxLine_2_ext[i + 1] = line_idx + 6;

        cmsh::model::geo::addLine(point_idx + 2, point_idx + 4, line_idx + 7);

        int Curve_tag1 = cmsh::model::geo::addCurveLoop({line_idx + 1, line_idx + 3, -(line_idx + 2), -LastyLine});
        cmsh::model::geo::addSurfaceFilling({Curve_tag1}, surface_idx + 1);

        int Curve_tagl = cmsh::model::geo::addCurveLoop({LastLefty, line_idx + 4, -(line_idx + 5), -(line_idx + 1)});
        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 2);

        int Curve_tagr = cmsh::model::geo::addCurveLoop({LastRighty, line_idx + 6, -(line_idx + 7), -(line_idx + 2)});
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 3);

        LastSurface.push_back(surface_idx + 1);
        LastLeftSurf.push_back(surface_idx + 2);
        LastRightSurf.push_back(surface_idx + 3);

        LastxPoint_1 = point_idx + 1, LastxPoint_1_ext = point_idx + 3;
        LastxPoint_2 = point_idx + 2, LastxPoint_2_ext = point_idx + 4;
        LastyLine = line_idx + 3, LastLefty = line_idx + 5, LastRighty = line_idx + 7;
        point_idx += 4;
        line_idx += 7;
        surface_idx += 3;
    }

    LeftXpoint[0][num_surf] = 2;
    RightXpoint[0][num_surf] = 4;

    LeftXpoint_ext[0][num_surf + 1] = 8;
    LeftXpoint_ext[0][num_surf + 2] = 10;
    RightXpoint_ext[0][num_surf + 1] = 14;
    RightXpoint_ext[0][num_surf + 2] = 16;

    cmsh::model::geo::addLine(LastxPoint_1, 2, line_idx + 1);
    LastxLine_1[num_surf - 1] = line_idx + 1;
    cmsh::model::geo::addLine(LastxPoint_2, 4, line_idx + 2);
    LastxLine_2[num_surf - 1] = line_idx + 2;

    int Curve_tag2 = cmsh::model::geo::addCurveLoop({line_idx + 1, 2, -(line_idx + 2), -LastyLine});
    cmsh::model::geo::addSurfaceFilling({Curve_tag2}, surface_idx + 1);

    LastSurface.push_back(surface_idx + 1);

    cmsh::model::geo::addLine(1, 6, line_idx + 3);
    cmsh::model::geo::addLine(3, 12, line_idx + 4);

    int Curve_tag3 = cmsh::model::geo::addCurveLoop({line_idx + 3, 3, -(line_idx + 4), -1});
    cmsh::model::geo::addSurfaceFilling({Curve_tag3}, surface_idx + 2);
    LastLine_x1 = line_idx + 3, LastLine_x2 = line_idx + 4;

    cmsh::model::geo::addLine(2, 9, line_idx + 5);
    cmsh::model::geo::addLine(4, 15, line_idx + 6);

    int Curve_tag4 = cmsh::model::geo::addCurveLoop({line_idx + 5, 4, -(line_idx + 6), -2});
    cmsh::model::geo::addSurfaceFilling({Curve_tag4}, surface_idx + 3);
    LastLine_x3 = line_idx + 5, LastLine_x4 = line_idx + 6;

    LastFrontSurf = surface_idx + 2;
    LastBackSurf = surface_idx + 3;

    cmsh::model::geo::addLine(LastxPoint_1_ext, 8, line_idx + 7);
    cmsh::model::geo::addLine(8, 10, line_idx + 8);
    LastxLine_1_ext[num_surf] = line_idx + 7;
    LastxLine_1_ext[num_surf + 1] = line_idx + 8;

    cmsh::model::geo::addLine(LastxPoint_2_ext, 14, line_idx + 9);
    cmsh::model::geo::addLine(14, 16, line_idx + 10);
    LastxLine_2_ext[num_surf] = line_idx + 9;
    LastxLine_2_ext[num_surf + 1] = line_idx + 10;

    cmsh::model::geo::addLine(2, 8, line_idx + 11);
    cmsh::model::geo::addLine(9, 10, line_idx + 12);
    LastLeftLine_2 = line_idx + 12;

    cmsh::model::geo::addLine(4, 14, line_idx + 13);
    cmsh::model::geo::addLine(15, 16, line_idx + 14);
    LastRightLine_2 = line_idx + 14;

    int Curve_tagl = cmsh::model::geo::addCurveLoop({LastLefty, line_idx + 7, -(line_idx + 11), -(line_idx + 1)});
    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 4);
    LastLeftSurf.push_back(surface_idx + 4);

    int Curve_tagr = cmsh::model::geo::addCurveLoop({LastRighty, line_idx + 9, -(line_idx + 13), -(line_idx + 2)});
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 5);
    LastRightSurf.push_back(surface_idx + 5);

    Curve_tagl = cmsh::model::geo::addCurveLoop({line_idx + 11, line_idx + 8, -(line_idx + 12), -(line_idx + 5)});
    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 6);
    LastLeftSurf.push_back(surface_idx + 6);

    Curve_tagr = cmsh::model::geo::addCurveLoop({line_idx + 13, line_idx + 10, -(line_idx + 14), -(line_idx + 6)});
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 7);
    LastRightSurf.push_back(surface_idx + 7);

    Curve_tagl = cmsh::model::geo::addCurveLoop({6, 5, -7, -(line_idx + 3)});
    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 8);
    LastLeftSurf.push_back(surface_idx + 8);

    Curve_tagr = cmsh::model::geo::addCurveLoop({9, 8, -10, -(line_idx + 4)});
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 9);
    LastRightSurf.push_back(surface_idx + 9);

    line_idx += 14;
    surface_idx += 9;

    LastPoint_1 = 1, LastPoint_2 = 2, LastPoint_3 = 3, LastPoint_4 = 4;
    LastPoint_5 = 5, LastPoint_6 = 6, LastPoint_7 = 7, LastPoint_8 = 8, LastPoint_9 = 9, LastPoint_10 = 10;
    LastPoint_11 = 11, LastPoint_12 = 12, LastPoint_13 = 13, LastPoint_14 = 14, LastPoint_15 = 15, LastPoint_16 = 16;
    LastLine_1 = 1, LastLine_2 = 2;

    EnterLastPoint.push_back(point_idx);

    Interface_surf[0].insert(Interface_surf[0].end(), LastSurface.begin(), LastSurface.end());
    Interface_surf[0].insert(Interface_surf[0].end(), LastLeftSurf.begin(), LastLeftSurf.end());
    Interface_surf[0].insert(Interface_surf[0].end(), LastRightSurf.begin(), LastRightSurf.end());
    Interface_surf[0].push_back(LastFrontSurf);
    Interface_surf[0].push_back(LastBackSurf);

    for (int i = 1; i <= line_idx; ++i) {
        Interface_line[0].push_back(i);
    }

    // 0 层
    int down = 0, up = 1;
    std::vector<int> surfaceTags;
    std::vector<int> FrontsurfaceTags, BacksurfaceTags;
    std::vector<int> LeftsurfaceTags, RightsurfaceTags;
    surfaceTags.assign(LastSurface.begin(), LastSurface.end());
    LeftsurfaceTags.assign(LastLeftSurf.begin(), LastLeftSurf.end());
    RightsurfaceTags.assign(LastRightSurf.begin(), LastRightSurf.end());
    FrontsurfaceTags.push_back(LastFrontSurf);
    BacksurfaceTags.push_back(LastBackSurf);

    LastSurface.clear();
    LastLeftSurf.clear();
    LastRightSurf.clear();

    compoundLines_y2.clear();
    compoundLines_y4.clear();

    bcIdx_1 = 0 * NX         + 0;
    bcIdx_2 = 0 * NX         + NX - 1;
    bcIdx_3 = (NY - 1) * NX + 0;
    bcIdx_4 = (NY - 1) * NX + NX - 1;

    int EndPoint_idx1 = LastPoint_1;
    int EndPoint_idx2 = point_idx + 1;
    int EndPoint_idx3 = LastPoint_2;
    int EndPoint_idx4 = point_idx + 2;
    int EndPoint_idx5 = LastPoint_3;
    int EndPoint_idx6 = point_idx + 3;
    int EndPoint_idx7 = LastPoint_4;
    int EndPoint_idx8 = point_idx + 4;

    int EndPoint_idx9 = point_idx + 5;
    int EndPoint_idx10 = point_idx + 6;
    int EndPoint_idx11 = point_idx + 7;

    int EndPoint_idx12 = point_idx + 8;
    int EndPoint_idx13 = point_idx + 9;
    int EndPoint_idx14 = point_idx + 10;

    int EndPoint_idx15 = point_idx + 11;
    int EndPoint_idx16 = point_idx + 12;
    int EndPoint_idx17 = point_idx + 13;

    int EndPoint_idx18 = point_idx + 14;
    int EndPoint_idx19 = point_idx + 15;
    int EndPoint_idx20 = point_idx + 16;

    int Line_idx1 = line_idx + 1;
    int Line_idx2 = line_idx + 2;
    int Line_idx3 = line_idx + 3;
    int Line_idx4 = line_idx + 4;
    int Line_idx5 = line_idx + 5;
    int Line_idx6 = line_idx + 6;

    int Line_idx7 = line_idx + 7;
    int Line_idx8 = line_idx + 8;
    int Line_idx9 = line_idx + 9;
    int Line_idx10 = line_idx + 10;
    int Line_idx11 = line_idx + 11;
    int Line_idx12 = line_idx + 12;
    int Line_idx13 = line_idx + 13;
    int Line_idx14 = line_idx + 14;
    int Line_idx15 = line_idx + 15;
    int Line_idx16 = line_idx + 16;
    int Line_idx17 = line_idx + 17;
    int Line_idx18 = line_idx + 18;

    int Line_idx19 = line_idx + 19;
    int Line_idx20 = line_idx + 20;

    int Line_idx21 = line_idx + 21;
    int Line_idx22 = line_idx + 22;
    int Line_idx23 = line_idx + 23;

    int Line_idx24 = line_idx + 24;
    int Line_idx25 = line_idx + 25;
    int Line_idx26 = line_idx + 26;

    int Line_idx27 = line_idx + 27;
    int Line_idx28 = line_idx + 28;
    int Line_idx29 = line_idx + 29;

    int Line_idx30 = line_idx + 30;
    int Line_idx31 = line_idx + 31;
    int Line_idx32 = line_idx + 32;

    cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[0]._x, CellArray[bcIdx_1]._points[0]._y, CellArray[bcIdx_1]._points[0]._z, lc, EndPoint_idx2);
    cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[1]._x, CellArray[bcIdx_2]._points[1]._y, CellArray[bcIdx_2]._points[1]._z, lc, EndPoint_idx4);
    cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[2]._x, CellArray[bcIdx_3]._points[2]._y, CellArray[bcIdx_3]._points[2]._z, lc, EndPoint_idx6);
    cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[3]._x, CellArray[bcIdx_4]._points[3]._y, CellArray[bcIdx_4]._points[3]._z, lc, EndPoint_idx8);

    cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[0]._x, Y_ext_max, CellArray[bcIdx_1]._points[0]._z, lc, EndPoint_idx9);
    cmsh::model::geo::addPoint(X_ext_min, CellArray[bcIdx_1]._points[0]._y, CellArray[bcIdx_1]._points[0]._z, lc, EndPoint_idx10);
    cmsh::model::geo::addPoint(X_ext_min, Y_ext_max, CellArray[bcIdx_1]._points[0]._z, lc, EndPoint_idx11);

    cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[1]._x, Y_ext_max, CellArray[bcIdx_2]._points[1]._z, lc, EndPoint_idx12);
    cmsh::model::geo::addPoint(X_ext_max, CellArray[bcIdx_2]._points[1]._y, CellArray[bcIdx_2]._points[1]._z, lc, EndPoint_idx13);
    cmsh::model::geo::addPoint(X_ext_max, Y_ext_max, CellArray[bcIdx_2]._points[1]._z, lc, EndPoint_idx14);

    cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[2]._x, Y_ext_min, CellArray[bcIdx_3]._points[2]._z, lc, EndPoint_idx15);
    cmsh::model::geo::addPoint(X_ext_min, CellArray[bcIdx_3]._points[2]._y, CellArray[bcIdx_3]._points[2]._z, lc, EndPoint_idx16);
    cmsh::model::geo::addPoint(X_ext_min, Y_ext_min, CellArray[bcIdx_3]._points[2]._z, lc, EndPoint_idx17);

    cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[3]._x, Y_ext_min, CellArray[bcIdx_4]._points[3]._z, lc, EndPoint_idx18);
    cmsh::model::geo::addPoint(X_ext_max, CellArray[bcIdx_4]._points[3]._y, CellArray[bcIdx_4]._points[3]._z, lc, EndPoint_idx19);
    cmsh::model::geo::addPoint(X_ext_max, Y_ext_min, CellArray[bcIdx_4]._points[3]._z, lc, EndPoint_idx20);

    compoundLines_z1.push_back(EndPoint_idx2);
    compoundLines_z2.push_back(EndPoint_idx4);
    compoundLines_z3.push_back(EndPoint_idx6);
    compoundLines_z4.push_back(EndPoint_idx8);
            
    cmsh::model::geo::addLine(EndPoint_idx1, EndPoint_idx2, Line_idx1);
    cmsh::model::geo::addLine(EndPoint_idx3, EndPoint_idx4, Line_idx2);
    cmsh::model::geo::addLine(EndPoint_idx5, EndPoint_idx6, Line_idx3);
    cmsh::model::geo::addLine(EndPoint_idx7, EndPoint_idx8, Line_idx4);

    cmsh::model::geo::addLine(LastPoint_5, EndPoint_idx9, Line_idx7);
    cmsh::model::geo::addLine(LastPoint_6, EndPoint_idx10, Line_idx8);
    cmsh::model::geo::addLine(LastPoint_7, EndPoint_idx11, Line_idx9);

    cmsh::model::geo::addLine(LastPoint_8, EndPoint_idx12, Line_idx10);
    cmsh::model::geo::addLine(LastPoint_9, EndPoint_idx13, Line_idx11);
    cmsh::model::geo::addLine(LastPoint_10, EndPoint_idx14, Line_idx12);

    cmsh::model::geo::addLine(LastPoint_11, EndPoint_idx15, Line_idx13);
    cmsh::model::geo::addLine(LastPoint_12, EndPoint_idx16, Line_idx14);
    cmsh::model::geo::addLine(LastPoint_13, EndPoint_idx17, Line_idx15);

    cmsh::model::geo::addLine(LastPoint_14, EndPoint_idx18, Line_idx16);
    cmsh::model::geo::addLine(LastPoint_15, EndPoint_idx19, Line_idx17);
    cmsh::model::geo::addLine(LastPoint_16, EndPoint_idx20, Line_idx18);

    point_idx += 16;

    compoundLines_y2.push_back(EndPoint_idx2);
    compoundLines_y4.push_back(EndPoint_idx4);

    LeftXpoint[up][0] = EndPoint_idx2;
    RightXpoint[up][0] = EndPoint_idx6;

    LeftXpoint_ext[up][0] = EndPoint_idx11;
    LeftXpoint_ext[up][1] = EndPoint_idx9;

    RightXpoint_ext[up][0] = EndPoint_idx17;
    RightXpoint_ext[up][1] = EndPoint_idx15;

    for (int i = 0; i < num_surf - 1; i ++) {
        // j = 0
        bcIdx_3 = 0 * NX + i;
        // j = NY - 1
        bcIdx_4 = (NY - 1) * NX + i;

        cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[1]._x, CellArray[bcIdx_3]._points[1]._y, CellArray[bcIdx_3]._points[1]._z, lc, point_idx + 1);
        LeftXpoint[up][i + 1] = point_idx + 1;
        cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[3]._x, CellArray[bcIdx_4]._points[3]._y, CellArray[bcIdx_4]._points[3]._z, lc, point_idx + 2);
        RightXpoint[up][i + 1] = point_idx + 2;
        cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[1]._x, Y_ext_max, CellArray[bcIdx_3]._points[1]._z, lc, point_idx + 3);
        LeftXpoint_ext[up][i + 2] = point_idx + 3;
        cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[3]._x, Y_ext_min, CellArray[bcIdx_4]._points[3]._z, lc, point_idx + 4);
        RightXpoint_ext[up][i + 2] = point_idx + 4;
        point_idx += 4;
    }

    for (int j = 0; j < NY - 1; j++) {
        // i = 0
        bcIdx_1 = j * NX + 0; 
        // i = NX - 1
        bcIdx_2 = j * NX + NX - 1;

        cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[2]._x, CellArray[bcIdx_1]._points[2]._y, CellArray[bcIdx_1]._points[2]._z, lc, point_idx + 1);
        compoundLines_y2.push_back(point_idx + 1);
        cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[3]._x, CellArray[bcIdx_2]._points[3]._y, CellArray[bcIdx_2]._points[3]._z, lc, point_idx + 2);
        compoundLines_y4.push_back(point_idx + 2);
        point_idx += 2;
    }

    compoundLines_y2.push_back(EndPoint_idx6);
    compoundLines_y4.push_back(EndPoint_idx8);
    LeftXpoint[up][num_surf] = EndPoint_idx4;
    RightXpoint[up][num_surf] = EndPoint_idx8;

    LeftXpoint_ext[up][num_surf + 1] = EndPoint_idx12;
    LeftXpoint_ext[up][num_surf + 2] = EndPoint_idx14;
    RightXpoint_ext[up][num_surf + 1] = EndPoint_idx18;
    RightXpoint_ext[up][num_surf + 2] = EndPoint_idx20;

    cmsh::model::geo::addBSpline(compoundLines_y2, Line_idx5);
    cmsh::model::geo::addBSpline(compoundLines_y4, Line_idx6);

    cmsh::model::geo::addLine(EndPoint_idx10, EndPoint_idx16, Line_idx19);
    cmsh::model::geo::addLine(EndPoint_idx13, EndPoint_idx19, Line_idx20);

    Interface_line[1].push_back(Line_idx5);
    Interface_line[1].push_back(Line_idx6);

    Interface_line[1].push_back(Line_idx19);
    Interface_line[1].push_back(Line_idx20);

    // point 1 (EndPoint_idx2)
    cmsh::model::geo::addLine(EndPoint_idx9, EndPoint_idx11, Line_idx21);

    cmsh::model::geo::addLine(EndPoint_idx10, EndPoint_idx11, Line_idx22);
    cmsh::model::geo::addLine(EndPoint_idx2, EndPoint_idx9, Line_idx23);

    // point 3 (EndPoint_idx6)
    cmsh::model::geo::addLine(EndPoint_idx15, EndPoint_idx17, Line_idx24);

    cmsh::model::geo::addLine(EndPoint_idx16, EndPoint_idx17, Line_idx25);
    cmsh::model::geo::addLine(EndPoint_idx6, EndPoint_idx15, Line_idx26);

    // point 2 (EndPoint_idx4)
    cmsh::model::geo::addLine(EndPoint_idx12, EndPoint_idx14, Line_idx27);

    cmsh::model::geo::addLine(EndPoint_idx13, EndPoint_idx14, Line_idx28);
    cmsh::model::geo::addLine(EndPoint_idx4, EndPoint_idx12, Line_idx29);

    // point 4 (EndPoint_idx8)
    cmsh::model::geo::addLine(EndPoint_idx18, EndPoint_idx20, Line_idx30);

    cmsh::model::geo::addLine(EndPoint_idx19, EndPoint_idx20, Line_idx31);
    cmsh::model::geo::addLine(EndPoint_idx8, EndPoint_idx18, Line_idx32);

    Interface_line[1].push_back(Line_idx21);
    Interface_line[1].push_back(Line_idx22);
    Interface_line[1].push_back(Line_idx23);

    Interface_line[1].push_back(Line_idx24);
    Interface_line[1].push_back(Line_idx25);
    Interface_line[1].push_back(Line_idx26);

    Interface_line[1].push_back(Line_idx27);
    Interface_line[1].push_back(Line_idx28);
    Interface_line[1].push_back(Line_idx29);

    Interface_line[1].push_back(Line_idx30);
    Interface_line[1].push_back(Line_idx31);
    Interface_line[1].push_back(Line_idx32);

    line_idx += 32;

    // addCurveLoop
    Curve_tag3 = cmsh::model::geo::addCurveLoop({Line_idx1, Line_idx5, -Line_idx3, -LastLine_1});
    Curve_tag4 = cmsh::model::geo::addCurveLoop({Line_idx2, Line_idx6, -Line_idx4, -LastLine_2});

    int Curve_tag_1 = cmsh::model::geo::addCurveLoop({Line_idx8, Line_idx19, -Line_idx14, -LastLine_3});
    int Curve_tag_2 = cmsh::model::geo::addCurveLoop({Line_idx11, Line_idx20, -Line_idx17, -LastLine_4});

    Curve_tagl = cmsh::model::geo::addCurveLoop({LastLeftLine_1, Line_idx9, -Line_idx22, -Line_idx8});
    Curve_tagr = cmsh::model::geo::addCurveLoop({LastRightLine_1, Line_idx15, -Line_idx25, -Line_idx14});

    // addPlaneSurface
    cmsh::model::geo::addSurfaceFilling({Curve_tag3}, surface_idx + 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tag4}, surface_idx + 2);
    cmsh::model::geo::addSurfaceFilling({Curve_tag_1}, surface_idx + 3);
    cmsh::model::geo::addSurfaceFilling({Curve_tag_2}, surface_idx + 4);
    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 5);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 6);

    Curve_tagl = cmsh::model::geo::addCurveLoop({LastLeftLine_2, Line_idx12, -Line_idx28, -Line_idx11});
    Curve_tagr = cmsh::model::geo::addCurveLoop({LastRightLine_2, Line_idx18, -Line_idx31, -Line_idx17});

    LastLeftLine_1 = Line_idx22, LastRightLine_1 = Line_idx25;
    LastLeftLine_2 = Line_idx28, LastRightLine_2 = Line_idx31;

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 7);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 8);

    surfaceTags.push_back(surface_idx + 1);
    surfaceTags.push_back(surface_idx + 2);

    FrontsurfaceTags.push_back(surface_idx + 1);
    FrontsurfaceTags.push_back(surface_idx + 3);
    BacksurfaceTags.push_back(surface_idx + 2);
    BacksurfaceTags.push_back(surface_idx + 4);
    LeftsurfaceTags.push_back(surface_idx + 5);
    LeftsurfaceTags.push_back(surface_idx + 7);

    RightsurfaceTags.push_back(surface_idx + 6);
    RightsurfaceTags.push_back(surface_idx + 8);

    surface_idx += 8;

    LastyLine = Line_idx5, LastLefty = Line_idx23, LastRighty = Line_idx26;
    int LastzLine_1 = Line_idx1, LastzLine_2 = Line_idx3, LastzLine_3 = Line_idx7, LastzLine_4 = Line_idx13;
    int LeftXsurf_base = surface_idx + 1, RightXsurf_base = surface_idx + num_surf + 1, Zsurf_base = surface_idx + num_surf * 2 + 1;
    int LeftXsurf_base_ext = surface_idx + num_surf * 3 + 1, RightXsurf_base_ext = surface_idx + num_surf * 4 + 1;
    int LeftZsurf_base_ext = surface_idx + num_surf * 5 + 1, RightZsurf_base_ext = surface_idx + num_surf * 6 + 1;
    for (int i = 0; i < num_surf - 1; i++) {
        cmsh::model::geo::addLine(LeftXpoint[up][i + 1], RightXpoint[up][i + 1], line_idx + 1);

        cmsh::model::geo::addLine(LeftXpoint[up][i], LeftXpoint[up][i + 1], line_idx + 2);
        cmsh::model::geo::addLine(RightXpoint[up][i], RightXpoint[up][i + 1], line_idx + 3);

        cmsh::model::geo::addLine(LeftXpoint[down][i + 1], LeftXpoint[up][i + 1], line_idx + 4);
        cmsh::model::geo::addLine(RightXpoint[down][i + 1], RightXpoint[up][i + 1], line_idx + 5);

        Interface_line[1].push_back(line_idx + 1);
        Interface_line[1].push_back(line_idx + 2);
        Interface_line[1].push_back(line_idx + 3);

        int Curve_tag5 = cmsh::model::geo::addCurveLoop({line_idx + 2, line_idx + 1, -(line_idx + 3), -LastyLine});
        int Curve_tag6 = cmsh::model::geo::addCurveLoop({LastzLine_1, line_idx + 2, -(line_idx + 4), -(LastxLine_1[i])});
        int Curve_tag7 = cmsh::model::geo::addCurveLoop({LastzLine_2, line_idx + 3, -(line_idx + 5), -(LastxLine_2[i])});

        LastxLine_1[i] = line_idx + 2;
        LastxLine_2[i] = line_idx + 3;

        // left_ext
        cmsh::model::geo::addLine(LeftXpoint[up][i + 1], LeftXpoint_ext[up][i + 2], line_idx + 6);

        cmsh::model::geo::addLine(LeftXpoint_ext[up][i + 1], LeftXpoint_ext[up][i + 2], line_idx + 7);
        cmsh::model::geo::addLine(LeftXpoint_ext[down][i + 2], LeftXpoint_ext[up][i + 2], line_idx + 8);
        // right_ext
        cmsh::model::geo::addLine(RightXpoint[up][i + 1], RightXpoint_ext[up][i + 2], line_idx + 9);

        cmsh::model::geo::addLine(RightXpoint_ext[up][i + 1], RightXpoint_ext[up][i + 2], line_idx + 10);
        cmsh::model::geo::addLine(RightXpoint_ext[down][i + 2], RightXpoint_ext[up][i + 2], line_idx + 11);

        Interface_line[1].push_back(line_idx + 6);
        Interface_line[1].push_back(line_idx + 7);

        Interface_line[1].push_back(line_idx + 9);
        Interface_line[1].push_back(line_idx + 10);

        Curve_tagl = cmsh::model::geo::addCurveLoop({LastLefty, line_idx + 7, -(line_idx + 6), -(line_idx + 2)});
        Curve_tagr = cmsh::model::geo::addCurveLoop({LastRighty, line_idx + 10, -(line_idx + 9), -(line_idx + 3)});

        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftZsurf_base_ext + i);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightZsurf_base_ext + i);

        Curve_tagl = cmsh::model::geo::addCurveLoop({LastzLine_3, line_idx + 7, -(line_idx + 8), -(LastxLine_1_ext[i + 1])});
        Curve_tagr = cmsh::model::geo::addCurveLoop({LastzLine_4, line_idx + 10, -(line_idx + 11), -(LastxLine_2_ext[i + 1])});

        LastxLine_1_ext[i + 1] = line_idx + 7;
        LastxLine_2_ext[i + 1] = line_idx + 10;

        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftXsurf_base_ext + i);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightXsurf_base_ext + i);

        cmsh::model::geo::addSurfaceFilling({Curve_tag5}, Zsurf_base + i);
        cmsh::model::geo::addSurfaceFilling({Curve_tag6}, LeftXsurf_base + i);
        cmsh::model::geo::addSurfaceFilling({Curve_tag7}, RightXsurf_base + i);

        surfaceTags.push_back(Zsurf_base + i);
        surfaceTags.push_back(LeftXsurf_base + i);
        surfaceTags.push_back(RightXsurf_base + i);

        LastSurface.push_back(Zsurf_base + i);
        LastLeftSurf.push_back(LeftZsurf_base_ext + i);
        LastRightSurf.push_back(RightZsurf_base_ext + i);

        LeftsurfaceTags.push_back(LeftXsurf_base + i);
        LeftsurfaceTags.push_back(LeftXsurf_base_ext + i);
        LeftsurfaceTags.push_back(LeftZsurf_base_ext + i);
        RightsurfaceTags.push_back(RightXsurf_base + i);
        RightsurfaceTags.push_back(RightXsurf_base_ext + i);
        RightsurfaceTags.push_back(RightZsurf_base_ext + i);

        LastyLine = line_idx + 1, LastLefty = line_idx + 6, LastRighty = line_idx + 9;
        LastzLine_1 = line_idx + 4;
        LastzLine_2 = line_idx + 5;
        LastzLine_3 = line_idx + 8;
        LastzLine_4 = line_idx + 11;
        line_idx += 11;
        surface_idx += 7;
    }

    cmsh::model::geo::addLine(LeftXpoint[up][num_surf - 1], LeftXpoint[up][num_surf], line_idx + 1);
    cmsh::model::geo::addLine(RightXpoint[up][num_surf - 1], RightXpoint[up][num_surf], line_idx + 2);

    Interface_line[1].push_back(line_idx + 1);
    Interface_line[1].push_back(line_idx + 2);

    int Curve_tag8 = cmsh::model::geo::addCurveLoop({line_idx + 1, Line_idx6, -(line_idx + 2), -LastyLine});
    int Curve_tag9 = cmsh::model::geo::addCurveLoop({LastzLine_1, line_idx + 1, -(Line_idx2), -(LastxLine_1[num_surf - 1])});
    int Curve_tag10 = cmsh::model::geo::addCurveLoop({LastzLine_2, line_idx + 2, -(Line_idx4), -(LastxLine_2[num_surf - 1])});

    LastxLine_1[num_surf - 1] = line_idx + 1;
    LastxLine_2[num_surf - 1] = line_idx + 2;

    cmsh::model::geo::addSurfaceFilling({Curve_tag8} , Zsurf_base + num_surf - 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tag9} , LeftXsurf_base + num_surf - 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tag10}, RightXsurf_base + num_surf - 1);

    surfaceTags.push_back(Zsurf_base + num_surf - 1);
    surfaceTags.push_back(LeftXsurf_base + num_surf - 1);
    surfaceTags.push_back(RightXsurf_base + num_surf - 1);

    LastSurface.push_back(Zsurf_base + num_surf - 1);

    LeftsurfaceTags.push_back(LeftXsurf_base + num_surf - 1);
    RightsurfaceTags.push_back(RightXsurf_base + num_surf - 1);

    cmsh::model::geo::addLine(EndPoint_idx2, EndPoint_idx10, line_idx + 3);
    cmsh::model::geo::addLine(EndPoint_idx6, EndPoint_idx16, line_idx + 4);

    cmsh::model::geo::addLine(EndPoint_idx4, EndPoint_idx13, line_idx + 5);
    cmsh::model::geo::addLine(EndPoint_idx8, EndPoint_idx19, line_idx + 6);

    // left_ext
    cmsh::model::geo::addLine(LeftXpoint_ext[up][num_surf], LeftXpoint_ext[up][num_surf + 1], line_idx + 7);

    // right_ext
    cmsh::model::geo::addLine(RightXpoint_ext[up][num_surf], RightXpoint_ext[up][num_surf + 1], line_idx + 8);

    Interface_line[1].push_back(line_idx + 3);
    Interface_line[1].push_back(line_idx + 4);
    Interface_line[1].push_back(line_idx + 5);
    Interface_line[1].push_back(line_idx + 6);
    Interface_line[1].push_back(line_idx + 7);
    Interface_line[1].push_back(line_idx + 8);

    Curve_tagl = cmsh::model::geo::addCurveLoop({LastLefty, line_idx + 7, -Line_idx29, -(line_idx + 1)});
    Curve_tagr = cmsh::model::geo::addCurveLoop({LastRighty, line_idx + 8, -Line_idx32, -(line_idx + 2)});

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftZsurf_base_ext + num_surf - 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightZsurf_base_ext + num_surf - 1);

    Curve_tagl = cmsh::model::geo::addCurveLoop({LastzLine_3, line_idx + 7, -Line_idx10, -(LastxLine_1_ext[num_surf])});
    Curve_tagr = cmsh::model::geo::addCurveLoop({LastzLine_4, line_idx + 8, -Line_idx16, -(LastxLine_2_ext[num_surf])});

    LastxLine_1_ext[num_surf] = line_idx + 7;
    LastxLine_2_ext[num_surf] = line_idx + 8;

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftXsurf_base_ext + num_surf - 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightXsurf_base_ext + num_surf - 1);

    LastLeftSurf.push_back(LeftZsurf_base_ext + num_surf - 1);
    LastRightSurf.push_back(RightZsurf_base_ext + num_surf - 1);

    LeftsurfaceTags.push_back(LeftXsurf_base_ext + num_surf - 1);
    LeftsurfaceTags.push_back(LeftZsurf_base_ext + num_surf - 1);

    RightsurfaceTags.push_back(RightXsurf_base_ext + num_surf - 1);
    RightsurfaceTags.push_back(RightZsurf_base_ext + num_surf - 1);

    // other surfaces
    surface_idx += 7;

    Curve_tag8 = cmsh::model::geo::addCurveLoop({line_idx + 3, Line_idx19, -(line_idx + 4), -Line_idx5});
    Curve_tag9 = cmsh::model::geo::addCurveLoop({Line_idx1, line_idx + 3, -Line_idx8, -LastLine_x1});
    Curve_tag10 = cmsh::model::geo::addCurveLoop({Line_idx3, line_idx + 4, -Line_idx14, -LastLine_x2});

    cmsh::model::geo::addSurfaceFilling({Curve_tag8}, surface_idx + 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tag9}, surface_idx + 2);
    cmsh::model::geo::addSurfaceFilling({Curve_tag10}, surface_idx + 3);

    FrontsurfaceTags.push_back(surface_idx + 1);
    FrontsurfaceTags.push_back(surface_idx + 2);
    FrontsurfaceTags.push_back(surface_idx + 3);

    LeftsurfaceTags.push_back(surface_idx + 2);
    RightsurfaceTags.push_back(surface_idx + 3);

    Curve_tag8 = cmsh::model::geo::addCurveLoop({line_idx + 5, Line_idx20, -(line_idx + 6), -Line_idx6});
    Curve_tag9 = cmsh::model::geo::addCurveLoop({Line_idx2, line_idx + 5, -Line_idx11, -LastLine_x3});
    Curve_tag10 = cmsh::model::geo::addCurveLoop({Line_idx4, line_idx + 6, -Line_idx17, -LastLine_x4});

    cmsh::model::geo::addSurfaceFilling({Curve_tag8}, surface_idx + 4);
    cmsh::model::geo::addSurfaceFilling({Curve_tag9}, surface_idx + 5);
    cmsh::model::geo::addSurfaceFilling({Curve_tag10}, surface_idx + 6);

    BacksurfaceTags.push_back(surface_idx + 4);
    BacksurfaceTags.push_back(surface_idx + 5);
    BacksurfaceTags.push_back(surface_idx + 6);

    LeftsurfaceTags.push_back(surface_idx + 5);
    RightsurfaceTags.push_back(surface_idx + 6);

    Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx29, Line_idx27, -Line_idx28, -(line_idx + 5)});
    Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx32, Line_idx30, -Line_idx31, -(line_idx + 6)});

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 7);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 8);

    Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx23, Line_idx21, -Line_idx22, -(line_idx + 3)});
    Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx26, Line_idx24, -Line_idx25, -(line_idx + 4)});

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 9);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 10);

    LeftsurfaceTags.push_back(surface_idx + 7);
    LeftsurfaceTags.push_back(surface_idx + 9);
    LastLeftSurf.push_back(surface_idx + 7);
    LastLeftSurf.push_back(surface_idx + 9);

    RightsurfaceTags.push_back(surface_idx + 8);
    RightsurfaceTags.push_back(surface_idx + 10);
    LastRightSurf.push_back(surface_idx + 8);
    LastRightSurf.push_back(surface_idx + 10);

    Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx10, Line_idx27, -Line_idx12, -(LastxLine_1_ext[num_surf + 1])});
    Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx16, Line_idx30, -Line_idx18, -(LastxLine_2_ext[num_surf + 1])});

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 11);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 12);

    Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx7, Line_idx21, -Line_idx9, -(LastxLine_1_ext[0])});
    Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx13, Line_idx24, -Line_idx15, -(LastxLine_2_ext[0])});

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 13);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 14);

    LeftsurfaceTags.push_back(surface_idx + 11);
    LeftsurfaceTags.push_back(surface_idx + 13);
    RightsurfaceTags.push_back(surface_idx + 12);
    RightsurfaceTags.push_back(surface_idx + 14);

    LastLine_x1 = line_idx + 3, LastLine_x2 = line_idx + 4;
    LastLine_x3 = line_idx + 5, LastLine_x4 = line_idx + 6;
    LastFrontSurf = surface_idx + 1;
    LastBackSurf = surface_idx + 4;

    LastxLine_1_ext[0] = Line_idx21;
    LastxLine_1_ext[num_surf + 1] = Line_idx27;
    LastxLine_2_ext[0] = Line_idx24;
    LastxLine_2_ext[num_surf + 1] = Line_idx30;

    line_idx += 8;
    surface_idx += 14;

    cmsh::model::geo::addSurfaceLoop(surfaceTags, volume_idx + 1);
    cmsh::model::geo::addVolume({volume_idx + 1}, volume_idx + 1);

    cmsh::model::geo::addSurfaceLoop(FrontsurfaceTags, volume_idx + 2);
    cmsh::model::geo::addVolume({volume_idx + 2}, volume_idx + 2);

    cmsh::model::geo::addSurfaceLoop(BacksurfaceTags, volume_idx + 3);
    cmsh::model::geo::addVolume({volume_idx + 3}, volume_idx + 3);

    cmsh::model::geo::addSurfaceLoop(LeftsurfaceTags, volume_idx + 4);
    cmsh::model::geo::addVolume({volume_idx + 4}, volume_idx + 4);

    cmsh::model::geo::addSurfaceLoop(RightsurfaceTags, volume_idx + 5);
    cmsh::model::geo::addVolume({volume_idx + 5}, volume_idx + 5);

    cellzone_id[4].push_back(volume_idx + 1);
    cellzone_id[5].push_back(volume_idx + 2);
    cellzone_id[5].push_back(volume_idx + 3);
    cellzone_id[5].push_back(volume_idx + 4);
    cellzone_id[5].push_back(volume_idx + 5);

    volume_idx += 5;
    LastPoint_1 = EndPoint_idx2, LastPoint_2 = EndPoint_idx4, LastPoint_3 = EndPoint_idx6, LastPoint_4 = EndPoint_idx8;
    LastPoint_5 = EndPoint_idx9, LastPoint_6 = EndPoint_idx10, LastPoint_7 = EndPoint_idx11, LastPoint_8 = EndPoint_idx12, LastPoint_9 = EndPoint_idx13, LastPoint_10 = EndPoint_idx14;
    LastPoint_11 = EndPoint_idx15, LastPoint_12 = EndPoint_idx16, LastPoint_13 = EndPoint_idx17, LastPoint_14 = EndPoint_idx18, LastPoint_15 = EndPoint_idx19, LastPoint_16 = EndPoint_idx20;
    LastLine_1 = Line_idx5, LastLine_2 = Line_idx6;
    LastLine_3 = Line_idx19, LastLine_4 = Line_idx20;

    EnterLastPoint.push_back(point_idx);

    Interface_surf[1].insert(Interface_surf[1].end(), LastSurface.begin(), LastSurface.end());
    Interface_surf[1].insert(Interface_surf[1].end(), LastLeftSurf.begin(), LastLeftSurf.end());
    Interface_surf[1].insert(Interface_surf[1].end(), LastRightSurf.begin(), LastRightSurf.end());
    Interface_surf[1].push_back(LastFrontSurf);
    Interface_surf[1].push_back(LastBackSurf);

    for (int k = 0; k < NZ; k++) {
        down = (k + 1) % 2, up = k % 2;
        // std::vector<int> surfaceTags;
        surfaceTags.clear();
        FrontsurfaceTags.clear();
        BacksurfaceTags.clear();
        LeftsurfaceTags.clear();
        RightsurfaceTags.clear();

        FrontsurfaceTags.push_back(LastFrontSurf);
        BacksurfaceTags.push_back(LastBackSurf);
        surfaceTags.assign(LastSurface.begin(), LastSurface.end());
        LeftsurfaceTags.assign(LastLeftSurf.begin(), LastLeftSurf.end());
        RightsurfaceTags.assign(LastRightSurf.begin(), LastRightSurf.end());

        LastSurface.clear();
        LastLeftSurf.clear();
        LastRightSurf.clear();

        compoundLines_y2.clear();
        compoundLines_y4.clear();

        std::vector<int> compoundLines_y2_ext, compoundLines_y4_ext;

        int LayerOffset = k * NY * NX;

        bcIdx_1 = LayerOffset + 0 * NX         + 0;
        bcIdx_2 = LayerOffset + 0 * NX         + NX - 1;
        bcIdx_3 = LayerOffset + (NY - 1) * NX + 0;
        bcIdx_4 = LayerOffset + (NY - 1) * NX + NX - 1;

        EndPoint_idx1 = LastPoint_1;
        EndPoint_idx2 = point_idx + 1;
        EndPoint_idx3 = LastPoint_2;
        EndPoint_idx4 = point_idx + 2;
        EndPoint_idx5 = LastPoint_3;
        EndPoint_idx6 = point_idx + 3;
        EndPoint_idx7 = LastPoint_4;
        EndPoint_idx8 = point_idx + 4;

        EndPoint_idx9 = point_idx + 5;
        EndPoint_idx10 = point_idx + 6;
        EndPoint_idx11 = point_idx + 7;

        EndPoint_idx12 = point_idx + 8;
        EndPoint_idx13 = point_idx + 9;
        EndPoint_idx14 = point_idx + 10;

        EndPoint_idx15 = point_idx + 11;
        EndPoint_idx16 = point_idx + 12;
        EndPoint_idx17 = point_idx + 13;

        EndPoint_idx18 = point_idx + 14;
        EndPoint_idx19 = point_idx + 15;
        EndPoint_idx20 = point_idx + 16;

        Line_idx1 = line_idx + 1;
        Line_idx2 = line_idx + 2;
        Line_idx3 = line_idx + 3;
        Line_idx4 = line_idx + 4;
        Line_idx5 = line_idx + 5;
        Line_idx6 = line_idx + 6;

        Line_idx7 = line_idx + 7;
        Line_idx8 = line_idx + 8;
        Line_idx9 = line_idx + 9;
        Line_idx10 = line_idx + 10;
        Line_idx11 = line_idx + 11;
        Line_idx12 = line_idx + 12;
        Line_idx13 = line_idx + 13;
        Line_idx14 = line_idx + 14;
        Line_idx15 = line_idx + 15;
        Line_idx16 = line_idx + 16;
        Line_idx17 = line_idx + 17;
        Line_idx18 = line_idx + 18;

        Line_idx19 = line_idx + 19;
        Line_idx20 = line_idx + 20;

        Line_idx21 = line_idx + 21;
        Line_idx22 = line_idx + 22;
        Line_idx23 = line_idx + 23;

        Line_idx24 = line_idx + 24;
        Line_idx25 = line_idx + 25;
        Line_idx26 = line_idx + 26;

        Line_idx27 = line_idx + 27;
        Line_idx28 = line_idx + 28;
        Line_idx29 = line_idx + 29;

        Line_idx30 = line_idx + 30;
        Line_idx31 = line_idx + 31;
        Line_idx32 = line_idx + 32;

        cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[4]._x, CellArray[bcIdx_1]._points[4]._y, CellArray[bcIdx_1]._points[4]._z, lc, EndPoint_idx2);
        cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[5]._x, CellArray[bcIdx_2]._points[5]._y, CellArray[bcIdx_2]._points[5]._z, lc, EndPoint_idx4);
        cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[6]._x, CellArray[bcIdx_3]._points[6]._y, CellArray[bcIdx_3]._points[6]._z, lc, EndPoint_idx6);
        cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[7]._x, CellArray[bcIdx_4]._points[7]._y, CellArray[bcIdx_4]._points[7]._z, lc, EndPoint_idx8);

        cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[4]._x, Y_ext_max, CellArray[bcIdx_1]._points[4]._z, lc, EndPoint_idx9);
        cmsh::model::geo::addPoint(X_ext_min, CellArray[bcIdx_1]._points[4]._y, CellArray[bcIdx_1]._points[4]._z, lc, EndPoint_idx10);
        cmsh::model::geo::addPoint(X_ext_min, Y_ext_max, CellArray[bcIdx_1]._points[4]._z, lc, EndPoint_idx11);

        cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[5]._x, Y_ext_max, CellArray[bcIdx_2]._points[5]._z, lc, EndPoint_idx12);
        cmsh::model::geo::addPoint(X_ext_max, CellArray[bcIdx_2]._points[5]._y, CellArray[bcIdx_2]._points[5]._z, lc, EndPoint_idx13);
        cmsh::model::geo::addPoint(X_ext_max, Y_ext_max, CellArray[bcIdx_2]._points[5]._z, lc, EndPoint_idx14);

        cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[6]._x, Y_ext_min, CellArray[bcIdx_3]._points[6]._z, lc, EndPoint_idx15);
        cmsh::model::geo::addPoint(X_ext_min, CellArray[bcIdx_3]._points[6]._y, CellArray[bcIdx_3]._points[6]._z, lc, EndPoint_idx16);
        cmsh::model::geo::addPoint(X_ext_min, Y_ext_min, CellArray[bcIdx_3]._points[6]._z, lc, EndPoint_idx17);

        cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[7]._x, Y_ext_min, CellArray[bcIdx_4]._points[7]._z, lc, EndPoint_idx18);
        cmsh::model::geo::addPoint(X_ext_max, CellArray[bcIdx_4]._points[7]._y, CellArray[bcIdx_4]._points[7]._z, lc, EndPoint_idx19);
        cmsh::model::geo::addPoint(X_ext_max, Y_ext_min, CellArray[bcIdx_4]._points[7]._z, lc, EndPoint_idx20);

        compoundLines_z1.push_back(EndPoint_idx2);
        compoundLines_z2.push_back(EndPoint_idx4);
        compoundLines_z3.push_back(EndPoint_idx6);
        compoundLines_z4.push_back(EndPoint_idx8);

        cmsh::model::geo::addLine(EndPoint_idx1, EndPoint_idx2, Line_idx1);
        cmsh::model::geo::addLine(EndPoint_idx3, EndPoint_idx4, Line_idx2);
        cmsh::model::geo::addLine(EndPoint_idx5, EndPoint_idx6, Line_idx3);
        cmsh::model::geo::addLine(EndPoint_idx7, EndPoint_idx8, Line_idx4);

        cmsh::model::geo::addLine(LastPoint_5, EndPoint_idx9, Line_idx7);
        cmsh::model::geo::addLine(LastPoint_6, EndPoint_idx10, Line_idx8);
        cmsh::model::geo::addLine(LastPoint_7, EndPoint_idx11, Line_idx9);

        cmsh::model::geo::addLine(LastPoint_8, EndPoint_idx12, Line_idx10);
        cmsh::model::geo::addLine(LastPoint_9, EndPoint_idx13, Line_idx11);
        cmsh::model::geo::addLine(LastPoint_10, EndPoint_idx14, Line_idx12);

        cmsh::model::geo::addLine(LastPoint_11, EndPoint_idx15, Line_idx13);
        cmsh::model::geo::addLine(LastPoint_12, EndPoint_idx16, Line_idx14);
        cmsh::model::geo::addLine(LastPoint_13, EndPoint_idx17, Line_idx15);

        cmsh::model::geo::addLine(LastPoint_14, EndPoint_idx18, Line_idx16);
        cmsh::model::geo::addLine(LastPoint_15, EndPoint_idx19, Line_idx17);
        cmsh::model::geo::addLine(LastPoint_16, EndPoint_idx20, Line_idx18);

        point_idx += 16;

        compoundLines_y2.push_back(EndPoint_idx2);
        compoundLines_y4.push_back(EndPoint_idx4);
        compoundLines_y2_ext.push_back(EndPoint_idx10);
        compoundLines_y4_ext.push_back(EndPoint_idx13);

        LeftXpoint[up][0] = EndPoint_idx2;
        RightXpoint[up][0] = EndPoint_idx6;

        LeftXpoint_ext[up][0] = EndPoint_idx11;
        LeftXpoint_ext[up][1] = EndPoint_idx9;

        RightXpoint_ext[up][0] = EndPoint_idx17;
        RightXpoint_ext[up][1] = EndPoint_idx15;

        for (int i = 0; i < num_surf - 1; i ++) {
            // j = 0
            bcIdx_3 = LayerOffset + 0 * NX + i;
            // j = NY - 1
            bcIdx_4 = LayerOffset + (NY - 1) * NX + i;

            cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[5]._x, CellArray[bcIdx_3]._points[5]._y, CellArray[bcIdx_3]._points[5]._z, lc, point_idx + 1);
            LeftXpoint[up][i + 1] = point_idx + 1;
            cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[7]._x, CellArray[bcIdx_4]._points[7]._y, CellArray[bcIdx_4]._points[7]._z, lc, point_idx + 2);
            RightXpoint[up][i + 1] = point_idx + 2;
            cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[5]._x, Y_ext_max, CellArray[bcIdx_3]._points[5]._z, lc, point_idx + 3);
            LeftXpoint_ext[up][i + 2] = point_idx + 3;
            cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[7]._x, Y_ext_min, CellArray[bcIdx_4]._points[7]._z, lc, point_idx + 4);
            RightXpoint_ext[up][i + 2] = point_idx + 4;
            point_idx += 4;
        }

        for (int j = 0; j < NY - 1; j++) {
            // i = 0
            bcIdx_1 = LayerOffset + j * NX + 0; 
            // i = NX - 1
            bcIdx_2 = LayerOffset + j * NX + NX - 1;

            cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[6]._x, CellArray[bcIdx_1]._points[6]._y, CellArray[bcIdx_1]._points[6]._z, lc, point_idx + 1);
            compoundLines_y2.push_back(point_idx + 1);
            cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[7]._x, CellArray[bcIdx_2]._points[7]._y, CellArray[bcIdx_2]._points[7]._z, lc, point_idx + 2);
            compoundLines_y4.push_back(point_idx + 2);
            cmsh::model::geo::addPoint(X_ext_min, CellArray[bcIdx_1]._points[6]._y, CellArray[bcIdx_1]._points[6]._z, lc, point_idx + 3);
            compoundLines_y2_ext.push_back(point_idx + 3);
            cmsh::model::geo::addPoint(X_ext_max, CellArray[bcIdx_2]._points[7]._y, CellArray[bcIdx_2]._points[7]._z, lc, point_idx + 4);
            compoundLines_y4_ext.push_back(point_idx + 4);
            point_idx += 4;
        }

        compoundLines_y2.push_back(EndPoint_idx6);
        compoundLines_y4.push_back(EndPoint_idx8);
        compoundLines_y2_ext.push_back(EndPoint_idx16);
        compoundLines_y4_ext.push_back(EndPoint_idx19);
        LeftXpoint[up][num_surf] = EndPoint_idx4;
        RightXpoint[up][num_surf] = EndPoint_idx8;

        LeftXpoint_ext[up][num_surf + 1] = EndPoint_idx12;
        LeftXpoint_ext[up][num_surf + 2] = EndPoint_idx14;
        RightXpoint_ext[up][num_surf + 1] = EndPoint_idx18;
        RightXpoint_ext[up][num_surf + 2] = EndPoint_idx20;

        cmsh::model::geo::addBSpline(compoundLines_y2, Line_idx5);
        cmsh::model::geo::addBSpline(compoundLines_y4, Line_idx6);

        cmsh::model::geo::addBSpline(compoundLines_y2_ext, Line_idx19);
        cmsh::model::geo::addBSpline(compoundLines_y4_ext, Line_idx20);

        Interface_line[k + 2].push_back(Line_idx5);
        Interface_line[k + 2].push_back(Line_idx6);

        Interface_line[k + 2].push_back(Line_idx19);
        Interface_line[k + 2].push_back(Line_idx20);

        // point 1 (EndPoint_idx2)
        cmsh::model::geo::addLine(EndPoint_idx9, EndPoint_idx11, Line_idx21);

        cmsh::model::geo::addLine(EndPoint_idx10, EndPoint_idx11, Line_idx22);
        cmsh::model::geo::addLine(EndPoint_idx2, EndPoint_idx9, Line_idx23);

        // point 3 (EndPoint_idx6)
        cmsh::model::geo::addLine(EndPoint_idx15, EndPoint_idx17, Line_idx24);

        cmsh::model::geo::addLine(EndPoint_idx16, EndPoint_idx17, Line_idx25);
        cmsh::model::geo::addLine(EndPoint_idx6, EndPoint_idx15, Line_idx26);

        // point 2 (EndPoint_idx4)
        cmsh::model::geo::addLine(EndPoint_idx12, EndPoint_idx14, Line_idx27);

        cmsh::model::geo::addLine(EndPoint_idx13, EndPoint_idx14, Line_idx28);
        cmsh::model::geo::addLine(EndPoint_idx4, EndPoint_idx12, Line_idx29);

        // point 4 (EndPoint_idx8)
        cmsh::model::geo::addLine(EndPoint_idx18, EndPoint_idx20, Line_idx30);

        cmsh::model::geo::addLine(EndPoint_idx19, EndPoint_idx20, Line_idx31);
        cmsh::model::geo::addLine(EndPoint_idx8, EndPoint_idx18, Line_idx32);

        Interface_line[k + 2].push_back(Line_idx21);
        Interface_line[k + 2].push_back(Line_idx22);
        Interface_line[k + 2].push_back(Line_idx23);

        Interface_line[k + 2].push_back(Line_idx24);
        Interface_line[k + 2].push_back(Line_idx25);
        Interface_line[k + 2].push_back(Line_idx26);

        Interface_line[k + 2].push_back(Line_idx27);
        Interface_line[k + 2].push_back(Line_idx28);
        Interface_line[k + 2].push_back(Line_idx29);

        Interface_line[k + 2].push_back(Line_idx30);
        Interface_line[k + 2].push_back(Line_idx31);
        Interface_line[k + 2].push_back(Line_idx32);

        line_idx += 32;

        // addCurveLoop
        Curve_tag3 = cmsh::model::geo::addCurveLoop({Line_idx1, Line_idx5, -Line_idx3, -LastLine_1});
        Curve_tag4 = cmsh::model::geo::addCurveLoop({Line_idx2, Line_idx6, -Line_idx4, -LastLine_2});

        Curve_tag_1 = cmsh::model::geo::addCurveLoop({Line_idx8, Line_idx19, -Line_idx14, -LastLine_3});
        Curve_tag_2 = cmsh::model::geo::addCurveLoop({Line_idx11, Line_idx20, -Line_idx17, -LastLine_4});

        Curve_tagl = cmsh::model::geo::addCurveLoop({LastLeftLine_1, Line_idx9, -Line_idx22, -Line_idx8});
        Curve_tagr = cmsh::model::geo::addCurveLoop({LastRightLine_1, Line_idx15, -Line_idx25, -Line_idx14});

        // addPlaneSurface
        cmsh::model::geo::addSurfaceFilling({Curve_tag3}, surface_idx + 1);
        cmsh::model::geo::addSurfaceFilling({Curve_tag4}, surface_idx + 2);
        cmsh::model::geo::addSurfaceFilling({Curve_tag_1}, surface_idx + 3);
        cmsh::model::geo::addSurfaceFilling({Curve_tag_2}, surface_idx + 4);
        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 5);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 6);

        Curve_tagl = cmsh::model::geo::addCurveLoop({LastLeftLine_2, Line_idx12, -Line_idx28, -Line_idx11});
        Curve_tagr = cmsh::model::geo::addCurveLoop({LastRightLine_2, Line_idx18, -Line_idx31, -Line_idx17});

        LastLeftLine_1 = Line_idx22, LastRightLine_1 = Line_idx25;
        LastLeftLine_2 = Line_idx28, LastRightLine_2 = Line_idx31;

        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 7);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 8);

        surfaceTags.push_back(surface_idx + 1);
        surfaceTags.push_back(surface_idx + 2);

        FrontsurfaceTags.push_back(surface_idx + 1);
        FrontsurfaceTags.push_back(surface_idx + 3);
        BacksurfaceTags.push_back(surface_idx + 2);
        BacksurfaceTags.push_back(surface_idx + 4);

        LeftsurfaceTags.push_back(surface_idx + 5);
        LeftsurfaceTags.push_back(surface_idx + 7);

        RightsurfaceTags.push_back(surface_idx + 6);
        RightsurfaceTags.push_back(surface_idx + 8);

        surface_idx += 8;

        LastyLine = Line_idx5, LastLefty = Line_idx23, LastRighty = Line_idx26;
        LastzLine_1 = Line_idx1, LastzLine_2 = Line_idx3, LastzLine_3 = Line_idx7, LastzLine_4 = Line_idx13;
        LeftXsurf_base = surface_idx + 1, RightXsurf_base = surface_idx + num_surf + 1, Zsurf_base = surface_idx + num_surf * 2 + 1;
        LeftXsurf_base_ext = surface_idx + num_surf * 3 + 1, RightXsurf_base_ext = surface_idx + num_surf * 4 + 1;
        LeftZsurf_base_ext = surface_idx + num_surf * 5 + 1, RightZsurf_base_ext = surface_idx + num_surf * 6 + 1;

        for (int i = 0; i < num_surf - 1; i++) {
            std::vector<int> tmp_compoundLines;
            tmp_compoundLines.resize(NY + 1);

            tmp_compoundLines[0] = LeftXpoint[up][i + 1];
            for (int j = 0; j < NY - 1; j++) {
                bcIdx_1 = LayerOffset + j * NX + i;
                cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[7]._x, CellArray[bcIdx_1]._points[7]._y, CellArray[bcIdx_1]._points[7]._z, lc, point_idx + 1);
                tmp_compoundLines[j + 1] = point_idx + 1;
                point_idx += 1; 
            }
            tmp_compoundLines[NY] = RightXpoint[up][i + 1];
            // printf("sizeof tmp_compoundLines: %d\n", tmp_compoundLines.size());
            cmsh::model::geo::addBSpline(tmp_compoundLines, line_idx + 1);
            // TODO
            cmsh::model::geo::addLine(LeftXpoint[up][i], LeftXpoint[up][i + 1], line_idx + 2);
            cmsh::model::geo::addLine(RightXpoint[up][i], RightXpoint[up][i + 1], line_idx + 3);

            cmsh::model::geo::addLine(LeftXpoint[down][i + 1], LeftXpoint[up][i + 1], line_idx + 4);
            cmsh::model::geo::addLine(RightXpoint[down][i + 1], RightXpoint[up][i + 1], line_idx + 5);

            Interface_line[k + 2].push_back(line_idx + 1);
            Interface_line[k + 2].push_back(line_idx + 2);
            Interface_line[k + 2].push_back(line_idx + 3);

            int Curve_tag5 = cmsh::model::geo::addCurveLoop({line_idx + 2, line_idx + 1, -(line_idx + 3), -LastyLine});
            int Curve_tag6 = cmsh::model::geo::addCurveLoop({LastzLine_1, line_idx + 2, -(line_idx + 4), -(LastxLine_1[i])});
            int Curve_tag7 = cmsh::model::geo::addCurveLoop({LastzLine_2, line_idx + 3, -(line_idx + 5), -(LastxLine_2[i])});

            LastxLine_1[i] = line_idx + 2;
            LastxLine_2[i] = line_idx + 3;

            // left_ext
            cmsh::model::geo::addLine(LeftXpoint[up][i + 1], LeftXpoint_ext[up][i + 2], line_idx + 6);

            cmsh::model::geo::addLine(LeftXpoint_ext[up][i + 1], LeftXpoint_ext[up][i + 2], line_idx + 7);
            cmsh::model::geo::addLine(LeftXpoint_ext[down][i + 2], LeftXpoint_ext[up][i + 2], line_idx + 8);
            // right_ext
            cmsh::model::geo::addLine(RightXpoint[up][i + 1], RightXpoint_ext[up][i + 2], line_idx + 9);

            cmsh::model::geo::addLine(RightXpoint_ext[up][i + 1], RightXpoint_ext[up][i + 2], line_idx + 10);
            cmsh::model::geo::addLine(RightXpoint_ext[down][i + 2], RightXpoint_ext[up][i + 2], line_idx + 11);

            Interface_line[k + 2].push_back(line_idx + 6);
            Interface_line[k + 2].push_back(line_idx + 7);

            Interface_line[k + 2].push_back(line_idx + 9);
            Interface_line[k + 2].push_back(line_idx + 10);

            Curve_tagl = cmsh::model::geo::addCurveLoop({LastLefty, line_idx + 7, -(line_idx + 6), -(line_idx + 2)});
            Curve_tagr = cmsh::model::geo::addCurveLoop({LastRighty, line_idx + 10, -(line_idx + 9), -(line_idx + 3)});

            cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftZsurf_base_ext + i);
            cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightZsurf_base_ext + i);

            Curve_tagl = cmsh::model::geo::addCurveLoop({LastzLine_3, line_idx + 7, -(line_idx + 8), -(LastxLine_1_ext[i + 1])});
            Curve_tagr = cmsh::model::geo::addCurveLoop({LastzLine_4, line_idx + 10, -(line_idx + 11), -(LastxLine_2_ext[i + 1])});

            LastxLine_1_ext[i + 1] = line_idx + 7;
            LastxLine_2_ext[i + 1] = line_idx + 10;

            cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftXsurf_base_ext + i);
            cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightXsurf_base_ext + i);

            cmsh::model::geo::addSurfaceFilling({Curve_tag5}, Zsurf_base + i);
            cmsh::model::geo::addSurfaceFilling({Curve_tag6}, LeftXsurf_base + i);
            cmsh::model::geo::addSurfaceFilling({Curve_tag7}, RightXsurf_base + i);

            surfaceTags.push_back(Zsurf_base + i);
            surfaceTags.push_back(LeftXsurf_base + i);
            surfaceTags.push_back(RightXsurf_base + i);

            LastSurface.push_back(Zsurf_base + i);
            LastLeftSurf.push_back(LeftZsurf_base_ext + i);
            LastRightSurf.push_back(RightZsurf_base_ext + i);

            LeftsurfaceTags.push_back(LeftXsurf_base + i);
            LeftsurfaceTags.push_back(LeftXsurf_base_ext + i);
            LeftsurfaceTags.push_back(LeftZsurf_base_ext + i);
            RightsurfaceTags.push_back(RightXsurf_base + i);
            RightsurfaceTags.push_back(RightXsurf_base_ext + i);
            RightsurfaceTags.push_back(RightZsurf_base_ext + i);

            LastyLine = line_idx + 1, LastLefty = line_idx + 6, LastRighty = line_idx + 9;
            LastzLine_1 = line_idx + 4;
            LastzLine_2 = line_idx + 5;
            LastzLine_3 = line_idx + 8;
            LastzLine_4 = line_idx + 11;
            line_idx += 11;
            surface_idx += 7;
        }

        cmsh::model::geo::addLine(LeftXpoint[up][num_surf - 1], LeftXpoint[up][num_surf], line_idx + 1);
        cmsh::model::geo::addLine(RightXpoint[up][num_surf - 1], RightXpoint[up][num_surf], line_idx + 2);

        Interface_line[k + 2].push_back(line_idx + 1);
        Interface_line[k + 2].push_back(line_idx + 2);

        Curve_tag8 = cmsh::model::geo::addCurveLoop({line_idx + 1, Line_idx6, -(line_idx + 2), -LastyLine});
        Curve_tag9 = cmsh::model::geo::addCurveLoop({LastzLine_1, line_idx + 1, -(Line_idx2), -(LastxLine_1[num_surf - 1])});
        Curve_tag10 = cmsh::model::geo::addCurveLoop({LastzLine_2, line_idx + 2, -(Line_idx4), -(LastxLine_2[num_surf - 1])});

        LastxLine_1[num_surf - 1] = line_idx + 1;
        LastxLine_2[num_surf - 1] = line_idx + 2;

        cmsh::model::geo::addSurfaceFilling({Curve_tag8} , Zsurf_base + num_surf - 1);
        cmsh::model::geo::addSurfaceFilling({Curve_tag9} , LeftXsurf_base + num_surf - 1);
        cmsh::model::geo::addSurfaceFilling({Curve_tag10}, RightXsurf_base + num_surf - 1);

        surfaceTags.push_back(Zsurf_base + num_surf - 1);
        surfaceTags.push_back(LeftXsurf_base + num_surf - 1);
        surfaceTags.push_back(RightXsurf_base + num_surf - 1);

        LastSurface.push_back(Zsurf_base + num_surf - 1);

        LeftsurfaceTags.push_back(LeftXsurf_base + num_surf - 1);
        RightsurfaceTags.push_back(RightXsurf_base + num_surf - 1);

        cmsh::model::geo::addLine(EndPoint_idx2, EndPoint_idx10, line_idx + 3);
        cmsh::model::geo::addLine(EndPoint_idx6, EndPoint_idx16, line_idx + 4);

        cmsh::model::geo::addLine(EndPoint_idx4, EndPoint_idx13, line_idx + 5);
        cmsh::model::geo::addLine(EndPoint_idx8, EndPoint_idx19, line_idx + 6);

        // left_ext
        cmsh::model::geo::addLine(LeftXpoint_ext[up][num_surf], LeftXpoint_ext[up][num_surf + 1], line_idx + 7);

        // right_ext
        cmsh::model::geo::addLine(RightXpoint_ext[up][num_surf], RightXpoint_ext[up][num_surf + 1], line_idx + 8);

        Interface_line[k + 2].push_back(line_idx + 3);
        Interface_line[k + 2].push_back(line_idx + 4);
        Interface_line[k + 2].push_back(line_idx + 5);
        Interface_line[k + 2].push_back(line_idx + 6);
        Interface_line[k + 2].push_back(line_idx + 7);
        Interface_line[k + 2].push_back(line_idx + 8);

        Curve_tagl = cmsh::model::geo::addCurveLoop({LastLefty, line_idx + 7, -Line_idx29, -(line_idx + 1)});
        Curve_tagr = cmsh::model::geo::addCurveLoop({LastRighty, line_idx + 8, -Line_idx32, -(line_idx + 2)});

        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftZsurf_base_ext + num_surf - 1);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightZsurf_base_ext + num_surf - 1);

        Curve_tagl = cmsh::model::geo::addCurveLoop({LastzLine_3, line_idx + 7, -Line_idx10, -(LastxLine_1_ext[num_surf])});
        Curve_tagr = cmsh::model::geo::addCurveLoop({LastzLine_4, line_idx + 8, -Line_idx16, -(LastxLine_2_ext[num_surf])});

        LastxLine_1_ext[num_surf] = line_idx + 7;
        LastxLine_2_ext[num_surf] = line_idx + 8;

        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftXsurf_base_ext + num_surf - 1);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightXsurf_base_ext + num_surf - 1);

        LastLeftSurf.push_back(LeftZsurf_base_ext + num_surf - 1);
        LastRightSurf.push_back(RightZsurf_base_ext + num_surf - 1);

        LeftsurfaceTags.push_back(LeftXsurf_base_ext + num_surf - 1);
        LeftsurfaceTags.push_back(LeftZsurf_base_ext + num_surf - 1);

        RightsurfaceTags.push_back(RightXsurf_base_ext + num_surf - 1);
        RightsurfaceTags.push_back(RightZsurf_base_ext + num_surf - 1);

        // other surfaces
        surface_idx += 7;

        Curve_tag8 = cmsh::model::geo::addCurveLoop({line_idx + 3, Line_idx19, -(line_idx + 4), -Line_idx5});
        Curve_tag9 = cmsh::model::geo::addCurveLoop({Line_idx1, line_idx + 3, -Line_idx8, -LastLine_x1});
        Curve_tag10 = cmsh::model::geo::addCurveLoop({Line_idx3, line_idx + 4, -Line_idx14, -LastLine_x2});

        cmsh::model::geo::addSurfaceFilling({Curve_tag8}, surface_idx + 1);
        cmsh::model::geo::addSurfaceFilling({Curve_tag9}, surface_idx + 2);
        cmsh::model::geo::addSurfaceFilling({Curve_tag10}, surface_idx + 3);

        FrontsurfaceTags.push_back(surface_idx + 1);
        FrontsurfaceTags.push_back(surface_idx + 2);
        FrontsurfaceTags.push_back(surface_idx + 3);

        LeftsurfaceTags.push_back(surface_idx + 2);
        RightsurfaceTags.push_back(surface_idx + 3);

        Curve_tag8 = cmsh::model::geo::addCurveLoop({line_idx + 5, Line_idx20, -(line_idx + 6), -Line_idx6});
        Curve_tag9 = cmsh::model::geo::addCurveLoop({Line_idx2, line_idx + 5, -Line_idx11, -LastLine_x3});
        Curve_tag10 = cmsh::model::geo::addCurveLoop({Line_idx4, line_idx + 6, -Line_idx17, -LastLine_x4});

        cmsh::model::geo::addSurfaceFilling({Curve_tag8}, surface_idx + 4);
        cmsh::model::geo::addSurfaceFilling({Curve_tag9}, surface_idx + 5);
        cmsh::model::geo::addSurfaceFilling({Curve_tag10}, surface_idx + 6);

        BacksurfaceTags.push_back(surface_idx + 4);
        BacksurfaceTags.push_back(surface_idx + 5);
        BacksurfaceTags.push_back(surface_idx + 6);

        LeftsurfaceTags.push_back(surface_idx + 5);
        RightsurfaceTags.push_back(surface_idx + 6);

        Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx29, Line_idx27, -Line_idx28, -(line_idx + 5)});
        Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx32, Line_idx30, -Line_idx31, -(line_idx + 6)});

        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 7);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 8);

        Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx23, Line_idx21, -Line_idx22, -(line_idx + 3)});
        Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx26, Line_idx24, -Line_idx25, -(line_idx + 4)});

        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 9);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 10);

        LeftsurfaceTags.push_back(surface_idx + 7);
        LeftsurfaceTags.push_back(surface_idx + 9);
        LastLeftSurf.push_back(surface_idx + 7);
        LastLeftSurf.push_back(surface_idx + 9);

        RightsurfaceTags.push_back(surface_idx + 8);
        RightsurfaceTags.push_back(surface_idx + 10);
        LastRightSurf.push_back(surface_idx + 8);
        LastRightSurf.push_back(surface_idx + 10);

        // std::cout << "Curve_tagr = " << Curve_tagr << std::endl;

        Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx10, Line_idx27, -Line_idx12, -(LastxLine_1_ext[num_surf + 1])});
        Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx16, Line_idx30, -Line_idx18, -(LastxLine_2_ext[num_surf + 1])});

        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 11);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 12);

        Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx7, Line_idx21, -Line_idx9, -(LastxLine_1_ext[0])});
        Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx13, Line_idx24, -Line_idx15, -(LastxLine_2_ext[0])});

        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 13);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 14);

        LeftsurfaceTags.push_back(surface_idx + 11);
        LeftsurfaceTags.push_back(surface_idx + 13);
        RightsurfaceTags.push_back(surface_idx + 12);
        RightsurfaceTags.push_back(surface_idx + 14);

        LastLine_x1 = line_idx + 3, LastLine_x2 = line_idx + 4;
        LastLine_x3 = line_idx + 5, LastLine_x4 = line_idx + 6;
        LastFrontSurf = surface_idx + 1;
        LastBackSurf = surface_idx + 4;

        LastxLine_1_ext[0] = Line_idx21;
        LastxLine_1_ext[num_surf + 1] = Line_idx27;
        LastxLine_2_ext[0] = Line_idx24;
        LastxLine_2_ext[num_surf + 1] = Line_idx30;

        line_idx += 8;
        surface_idx += 14;

        cmsh::model::geo::addSurfaceLoop(surfaceTags, volume_idx + 1);
        cmsh::model::geo::addVolume({volume_idx + 1}, volume_idx + 1);

        cmsh::model::geo::addSurfaceLoop(FrontsurfaceTags, volume_idx + 2);
        cmsh::model::geo::addVolume({volume_idx + 2}, volume_idx + 2);

        cmsh::model::geo::addSurfaceLoop(BacksurfaceTags, volume_idx + 3);
        cmsh::model::geo::addVolume({volume_idx + 3}, volume_idx + 3);

        cmsh::model::geo::addSurfaceLoop(LeftsurfaceTags, volume_idx + 4);
        cmsh::model::geo::addVolume({volume_idx + 4}, volume_idx + 4);

        cmsh::model::geo::addSurfaceLoop(RightsurfaceTags, volume_idx + 5);
        cmsh::model::geo::addVolume({volume_idx + 5}, volume_idx + 5);

        cellzone_id[0].push_back(volume_idx + 1);
        cellzone_id[1].push_back(volume_idx + 2);
        cellzone_id[1].push_back(volume_idx + 3);
        cellzone_id[1].push_back(volume_idx + 4);
        cellzone_id[1].push_back(volume_idx + 5);

        volume_idx += 5;
        LastPoint_1 = EndPoint_idx2, LastPoint_2 = EndPoint_idx4, LastPoint_3 = EndPoint_idx6, LastPoint_4 = EndPoint_idx8;
        LastPoint_5 = EndPoint_idx9, LastPoint_6 = EndPoint_idx10, LastPoint_7 = EndPoint_idx11, LastPoint_8 = EndPoint_idx12, LastPoint_9 = EndPoint_idx13, LastPoint_10 = EndPoint_idx14;
        LastPoint_11 = EndPoint_idx15, LastPoint_12 = EndPoint_idx16, LastPoint_13 = EndPoint_idx17, LastPoint_14 = EndPoint_idx18, LastPoint_15 = EndPoint_idx19, LastPoint_16 = EndPoint_idx20;
        LastLine_1 = Line_idx5, LastLine_2 = Line_idx6;
        LastLine_3 = Line_idx19, LastLine_4 = Line_idx20;

        EnterLastPoint.push_back(point_idx);

        Interface_surf[k + 2].insert(Interface_surf[k + 2].end(), LastSurface.begin(), LastSurface.end());
        Interface_surf[k + 2].insert(Interface_surf[k + 2].end(), LastLeftSurf.begin(), LastLeftSurf.end());
        Interface_surf[k + 2].insert(Interface_surf[k + 2].end(), LastRightSurf.begin(), LastRightSurf.end());
        Interface_surf[k + 2].push_back(LastFrontSurf);
        Interface_surf[k + 2].push_back(LastBackSurf);
    }

    cmsh::model::geo::synchronize();

    /* 2、Create Fault */ 
    std::vector<int> FaultSurface_idx;
    FaultSurface_idx.resize(FaultArray.size());

    for (int i = 0; i < FaultArray.size(); i++) {
        auto fault = FaultArray[i];
        std::vector<int> FaultCompoundLines_z1, FaultCompoundLines_z2;
        std::vector<int> FaultCompoundLines_l1, FaultCompoundLines_l2;
        std::vector<int> controlPoints;

        auto line = fault._lines[0];
        for (int j = 0; j < line.size(); j++) {
            auto point = line[j];
            cmsh::model::geo::addPoint(point._x, point._y, point._z, lc, point_idx + 1);
            FaultCompoundLines_z1.push_back(point_idx + 1);
            if (j == 0) {
                FaultCompoundLines_l1.push_back(point_idx + 1);
            } else if (j == line.size() - 1) {
                FaultCompoundLines_l2.push_back(point_idx + 1);
            }
            point_idx += 1;
        }

        for (int j = 1 ; j < fault._lines.size() - 1; j++) {
            auto line = fault._lines[j];
            for (int k = 0; k < line.size(); k++) {
                auto point = line[k];
                cmsh::model::geo::addPoint(point._x, point._y, point._z, lc, point_idx + 1);
                if (k == 0) {
                    FaultCompoundLines_l1.push_back(point_idx + 1);
                } else if (k == line.size() - 1) {
                    FaultCompoundLines_l2.push_back(point_idx + 1);
                } else {
                    controlPoints.push_back(point_idx + 1);
                }
                point_idx += 1;
            }
        }

        line = fault._lines[fault._lines.size() - 1];
        for (int j = 0; j < line.size(); j++) {
            auto point = line[j];
            cmsh::model::geo::addPoint(point._x, point._y, point._z, lc, point_idx + 1);
            FaultCompoundLines_z2.push_back(point_idx + 1);
            if (j == 0) {
                FaultCompoundLines_l1.push_back(point_idx + 1);
            } else if (j == line.size() - 1) {
                FaultCompoundLines_l2.push_back(point_idx + 1);
            }
            point_idx += 1;
        }

        cmsh::model::geo::addSpline(FaultCompoundLines_z1, line_idx + 1);
        cmsh::model::geo::addSpline(FaultCompoundLines_z2, line_idx + 2);
        cmsh::model::geo::addSpline(FaultCompoundLines_l1, line_idx + 3);
        cmsh::model::geo::addSpline(FaultCompoundLines_l2, line_idx + 4);

        cmsh::model::geo::addCurveLoop({line_idx + 1, line_idx + 4, -(line_idx + 2), -(line_idx + 3)}, surface_idx + 1);
        cmsh::model::geo::addSurfaceFilling({surface_idx + 1}, surface_idx + 1);
        FaultSurface_idx[i] = surface_idx + 1;

        surface_idx += 1;
        line_idx += 4;
    }

    /* 3、Grid extension */ 
    // Oversides
    int k = NZ;
    down = (k + 1) % 2, up = k % 2;
    surfaceTags.clear();
    FrontsurfaceTags.clear();
    BacksurfaceTags.clear();
    LeftsurfaceTags.clear();
    RightsurfaceTags.clear();

    FrontsurfaceTags.push_back(LastFrontSurf);
    BacksurfaceTags.push_back(LastBackSurf);
    surfaceTags.assign(LastSurface.begin(), LastSurface.end());
    LeftsurfaceTags.assign(LastLeftSurf.begin(), LastLeftSurf.end());
    RightsurfaceTags.assign(LastRightSurf.begin(), LastRightSurf.end());

    LastSurface.clear();
    LastLeftSurf.clear();
    LastRightSurf.clear();

    compoundLines_y2.clear();
    compoundLines_y4.clear();

    std::vector<int> compoundLines_y2_ext, compoundLines_y4_ext;

    int LayerOffset = (k - 1) * NY * NX;

    bcIdx_1 = LayerOffset + 0 * NX         + 0;
    bcIdx_2 = LayerOffset + 0 * NX         + NX - 1;
    bcIdx_3 = LayerOffset + (NY - 1) * NX + 0;
    bcIdx_4 = LayerOffset + (NY - 1) * NX + NX - 1;

    EndPoint_idx1 = LastPoint_1;
    EndPoint_idx2 = point_idx + 1;
    EndPoint_idx3 = LastPoint_2;
    EndPoint_idx4 = point_idx + 2;
    EndPoint_idx5 = LastPoint_3;
    EndPoint_idx6 = point_idx + 3;
    EndPoint_idx7 = LastPoint_4;
    EndPoint_idx8 = point_idx + 4;

    EndPoint_idx9 = point_idx + 5;
    EndPoint_idx10 = point_idx + 6;
    EndPoint_idx11 = point_idx + 7;

    EndPoint_idx12 = point_idx + 8;
    EndPoint_idx13 = point_idx + 9;
    EndPoint_idx14 = point_idx + 10;

    EndPoint_idx15 = point_idx + 11;
    EndPoint_idx16 = point_idx + 12;
    EndPoint_idx17 = point_idx + 13;

    EndPoint_idx18 = point_idx + 14;
    EndPoint_idx19 = point_idx + 15;
    EndPoint_idx20 = point_idx + 16;

    Line_idx1 = line_idx + 1;
    Line_idx2 = line_idx + 2;
    Line_idx3 = line_idx + 3;
    Line_idx4 = line_idx + 4;
    Line_idx5 = line_idx + 5;
    Line_idx6 = line_idx + 6;

    Line_idx7 = line_idx + 7;
    Line_idx8 = line_idx + 8;
    Line_idx9 = line_idx + 9;
    Line_idx10 = line_idx + 10;
    Line_idx11 = line_idx + 11;
    Line_idx12 = line_idx + 12;
    Line_idx13 = line_idx + 13;
    Line_idx14 = line_idx + 14;
    Line_idx15 = line_idx + 15;
    Line_idx16 = line_idx + 16;
    Line_idx17 = line_idx + 17;
    Line_idx18 = line_idx + 18;

    Line_idx19 = line_idx + 19;
    Line_idx20 = line_idx + 20;

    Line_idx21 = line_idx + 21;
    Line_idx22 = line_idx + 22;
    Line_idx23 = line_idx + 23;

    Line_idx24 = line_idx + 24;
    Line_idx25 = line_idx + 25;
    Line_idx26 = line_idx + 26;

    Line_idx27 = line_idx + 27;
    Line_idx28 = line_idx + 28;
    Line_idx29 = line_idx + 29;

    Line_idx30 = line_idx + 30;
    Line_idx31 = line_idx + 31;
    Line_idx32 = line_idx + 32;

    cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[4]._x, CellArray[bcIdx_1]._points[4]._y, Z_ext_max, lc, EndPoint_idx2);
    cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[5]._x, CellArray[bcIdx_2]._points[5]._y, Z_ext_max, lc, EndPoint_idx4);
    cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[6]._x, CellArray[bcIdx_3]._points[6]._y, Z_ext_max, lc, EndPoint_idx6);
    cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[7]._x, CellArray[bcIdx_4]._points[7]._y, Z_ext_max, lc, EndPoint_idx8);

    cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[4]._x, Y_ext_max, Z_ext_max, lc, EndPoint_idx9);
    cmsh::model::geo::addPoint(X_ext_min, CellArray[bcIdx_1]._points[4]._y, Z_ext_max, lc, EndPoint_idx10);
    cmsh::model::geo::addPoint(X_ext_min, Y_ext_max, Z_ext_max, lc, EndPoint_idx11);

    cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[5]._x, Y_ext_max, Z_ext_max, lc, EndPoint_idx12);
    cmsh::model::geo::addPoint(X_ext_max, CellArray[bcIdx_2]._points[5]._y, Z_ext_max, lc, EndPoint_idx13);
    cmsh::model::geo::addPoint(X_ext_max, Y_ext_max, Z_ext_max, lc, EndPoint_idx14);

    cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[6]._x, Y_ext_min, Z_ext_max, lc, EndPoint_idx15);
    cmsh::model::geo::addPoint(X_ext_min, CellArray[bcIdx_3]._points[6]._y, Z_ext_max, lc, EndPoint_idx16);
    cmsh::model::geo::addPoint(X_ext_min, Y_ext_min, Z_ext_max, lc, EndPoint_idx17);

    cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[7]._x, Y_ext_min, Z_ext_max, lc, EndPoint_idx18);
    cmsh::model::geo::addPoint(X_ext_max, CellArray[bcIdx_4]._points[7]._y, Z_ext_max, lc, EndPoint_idx19);
    cmsh::model::geo::addPoint(X_ext_max, Y_ext_min, Z_ext_max, lc, EndPoint_idx20);

    compoundLines_z1.push_back(EndPoint_idx2);
    compoundLines_z2.push_back(EndPoint_idx4);
    compoundLines_z3.push_back(EndPoint_idx6);
    compoundLines_z4.push_back(EndPoint_idx8);

    cmsh::model::geo::addLine(EndPoint_idx1, EndPoint_idx2, Line_idx1);
    cmsh::model::geo::addLine(EndPoint_idx3, EndPoint_idx4, Line_idx2);
    cmsh::model::geo::addLine(EndPoint_idx5, EndPoint_idx6, Line_idx3);
    cmsh::model::geo::addLine(EndPoint_idx7, EndPoint_idx8, Line_idx4);

    cmsh::model::geo::addLine(LastPoint_5, EndPoint_idx9, Line_idx7);
    cmsh::model::geo::addLine(LastPoint_6, EndPoint_idx10, Line_idx8);
    cmsh::model::geo::addLine(LastPoint_7, EndPoint_idx11, Line_idx9);

    cmsh::model::geo::addLine(LastPoint_8, EndPoint_idx12, Line_idx10);
    cmsh::model::geo::addLine(LastPoint_9, EndPoint_idx13, Line_idx11);
    cmsh::model::geo::addLine(LastPoint_10, EndPoint_idx14, Line_idx12);

    cmsh::model::geo::addLine(LastPoint_11, EndPoint_idx15, Line_idx13);
    cmsh::model::geo::addLine(LastPoint_12, EndPoint_idx16, Line_idx14);
    cmsh::model::geo::addLine(LastPoint_13, EndPoint_idx17, Line_idx15);

    cmsh::model::geo::addLine(LastPoint_14, EndPoint_idx18, Line_idx16);
    cmsh::model::geo::addLine(LastPoint_15, EndPoint_idx19, Line_idx17);
    cmsh::model::geo::addLine(LastPoint_16, EndPoint_idx20, Line_idx18);

    point_idx += 16;

    compoundLines_y2.push_back(EndPoint_idx2);
    compoundLines_y4.push_back(EndPoint_idx4);
    compoundLines_y2_ext.push_back(EndPoint_idx10);
    compoundLines_y4_ext.push_back(EndPoint_idx13);

    LeftXpoint[up][0] = EndPoint_idx2;
    RightXpoint[up][0] = EndPoint_idx6;

    LeftXpoint_ext[up][0] = EndPoint_idx11;
    LeftXpoint_ext[up][1] = EndPoint_idx9;

    RightXpoint_ext[up][0] = EndPoint_idx17;
    RightXpoint_ext[up][1] = EndPoint_idx15;

    for (int i = 0; i < num_surf - 1; i ++) {
        // j = 0
        bcIdx_3 = LayerOffset + 0 * NX + i;
        // j = NY - 1
        bcIdx_4 = LayerOffset + (NY - 1) * NX + i;

        cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[5]._x, CellArray[bcIdx_3]._points[5]._y, Z_ext_max, lc, point_idx + 1);
        LeftXpoint[up][i + 1] = point_idx + 1;
        cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[7]._x, CellArray[bcIdx_4]._points[7]._y, Z_ext_max, lc, point_idx + 2);
        RightXpoint[up][i + 1] = point_idx + 2;
        cmsh::model::geo::addPoint(CellArray[bcIdx_3]._points[5]._x, Y_ext_max, Z_ext_max, lc, point_idx + 3);
        LeftXpoint_ext[up][i + 2] = point_idx + 3;
        cmsh::model::geo::addPoint(CellArray[bcIdx_4]._points[7]._x, Y_ext_min, Z_ext_max, lc, point_idx + 4);
        RightXpoint_ext[up][i + 2] = point_idx + 4;
        point_idx += 4;
    }

    for (int j = 0; j < NY - 1; j++) {
        // i = 0
        bcIdx_1 = LayerOffset + j * NX + 0; 
        // i = NX - 1
        bcIdx_2 = LayerOffset + j * NX + NX - 1;

        cmsh::model::geo::addPoint(CellArray[bcIdx_1]._points[6]._x, CellArray[bcIdx_1]._points[6]._y, Z_ext_max, lc, point_idx + 1);
        compoundLines_y2.push_back(point_idx + 1);
        cmsh::model::geo::addPoint(CellArray[bcIdx_2]._points[7]._x, CellArray[bcIdx_2]._points[7]._y, Z_ext_max, lc, point_idx + 2);
        compoundLines_y4.push_back(point_idx + 2);
        point_idx += 2;
    }

    compoundLines_y2.push_back(EndPoint_idx6);
    compoundLines_y4.push_back(EndPoint_idx8);
    LeftXpoint[up][num_surf] = EndPoint_idx4;
    RightXpoint[up][num_surf] = EndPoint_idx8;

    LeftXpoint_ext[up][num_surf + 1] = EndPoint_idx12;
    LeftXpoint_ext[up][num_surf + 2] = EndPoint_idx14;
    RightXpoint_ext[up][num_surf + 1] = EndPoint_idx18;
    RightXpoint_ext[up][num_surf + 2] = EndPoint_idx20;

    cmsh::model::geo::addBSpline(compoundLines_y2, Line_idx5);
    cmsh::model::geo::addBSpline(compoundLines_y4, Line_idx6);

    cmsh::model::geo::addLine(EndPoint_idx10, EndPoint_idx16, Line_idx19);
    cmsh::model::geo::addLine(EndPoint_idx13, EndPoint_idx19, Line_idx20);

    Interface_line[NZ + 2].push_back(Line_idx5);
    Interface_line[NZ + 2].push_back(Line_idx6);

    Interface_line[NZ + 2].push_back(Line_idx19);
    Interface_line[NZ + 2].push_back(Line_idx20);

    // point 1 (EndPoint_idx2)
    cmsh::model::geo::addLine(EndPoint_idx9, EndPoint_idx11, Line_idx21);

    cmsh::model::geo::addLine(EndPoint_idx10, EndPoint_idx11, Line_idx22);
    cmsh::model::geo::addLine(EndPoint_idx2, EndPoint_idx9, Line_idx23);

    // point 3 (EndPoint_idx6)
    cmsh::model::geo::addLine(EndPoint_idx15, EndPoint_idx17, Line_idx24);

    cmsh::model::geo::addLine(EndPoint_idx16, EndPoint_idx17, Line_idx25);
    cmsh::model::geo::addLine(EndPoint_idx6, EndPoint_idx15, Line_idx26);

    // point 2 (EndPoint_idx4)
    cmsh::model::geo::addLine(EndPoint_idx12, EndPoint_idx14, Line_idx27);

    cmsh::model::geo::addLine(EndPoint_idx13, EndPoint_idx14, Line_idx28);
    cmsh::model::geo::addLine(EndPoint_idx4, EndPoint_idx12, Line_idx29);

    // point 4 (EndPoint_idx8)
    cmsh::model::geo::addLine(EndPoint_idx18, EndPoint_idx20, Line_idx30);

    cmsh::model::geo::addLine(EndPoint_idx19, EndPoint_idx20, Line_idx31);
    cmsh::model::geo::addLine(EndPoint_idx8, EndPoint_idx18, Line_idx32);

    Interface_line[NZ + 2].push_back(Line_idx21);
    Interface_line[NZ + 2].push_back(Line_idx22);
    Interface_line[NZ + 2].push_back(Line_idx23);

    Interface_line[NZ + 2].push_back(Line_idx24);
    Interface_line[NZ + 2].push_back(Line_idx25);
    Interface_line[NZ + 2].push_back(Line_idx26);

    Interface_line[NZ + 2].push_back(Line_idx27);
    Interface_line[NZ + 2].push_back(Line_idx28);
    Interface_line[NZ + 2].push_back(Line_idx29);

    Interface_line[NZ + 2].push_back(Line_idx30);
    Interface_line[NZ + 2].push_back(Line_idx31);
    Interface_line[NZ + 2].push_back(Line_idx32);

    line_idx += 32;

    // addCurveLoop
    Curve_tag3 = cmsh::model::geo::addCurveLoop({Line_idx1, Line_idx5, -Line_idx3, -LastLine_1});
    Curve_tag4 = cmsh::model::geo::addCurveLoop({Line_idx2, Line_idx6, -Line_idx4, -LastLine_2});

    Curve_tag_1 = cmsh::model::geo::addCurveLoop({Line_idx8, Line_idx19, -Line_idx14, -LastLine_3});
    Curve_tag_2 = cmsh::model::geo::addCurveLoop({Line_idx11, Line_idx20, -Line_idx17, -LastLine_4});

    Curve_tagl = cmsh::model::geo::addCurveLoop({LastLeftLine_1, Line_idx9, -Line_idx22, -Line_idx8});
    Curve_tagr = cmsh::model::geo::addCurveLoop({LastRightLine_1, Line_idx15, -Line_idx25, -Line_idx14});

    // addPlaneSurface
    cmsh::model::geo::addSurfaceFilling({Curve_tag3}, surface_idx + 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tag4}, surface_idx + 2);
    cmsh::model::geo::addSurfaceFilling({Curve_tag_1}, surface_idx + 3);
    cmsh::model::geo::addSurfaceFilling({Curve_tag_2}, surface_idx + 4);
    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 5);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 6);

    Curve_tagl = cmsh::model::geo::addCurveLoop({LastLeftLine_2, Line_idx12, -Line_idx28, -Line_idx11});
    Curve_tagr = cmsh::model::geo::addCurveLoop({LastRightLine_2, Line_idx18, -Line_idx31, -Line_idx17});

    LastLeftLine_1 = Line_idx22, LastRightLine_1 = Line_idx25;
    LastLeftLine_2 = Line_idx28, LastRightLine_2 = Line_idx31;

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 7);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 8);

    surfaceTags.push_back(surface_idx + 1);
    surfaceTags.push_back(surface_idx + 2);

    FrontsurfaceTags.push_back(surface_idx + 1);
    FrontsurfaceTags.push_back(surface_idx + 3);
    BacksurfaceTags.push_back(surface_idx + 2);
    BacksurfaceTags.push_back(surface_idx + 4);

    LeftsurfaceTags.push_back(surface_idx + 5);
    LeftsurfaceTags.push_back(surface_idx + 7);

    RightsurfaceTags.push_back(surface_idx + 6);
    RightsurfaceTags.push_back(surface_idx + 8);

    surface_idx += 8;

    LastyLine = Line_idx5, LastLefty = Line_idx23, LastRighty = Line_idx26;
    LastzLine_1 = Line_idx1, LastzLine_2 = Line_idx3, LastzLine_3 = Line_idx7, LastzLine_4 = Line_idx13;
    LeftXsurf_base = surface_idx + 1, RightXsurf_base = surface_idx + num_surf + 1, Zsurf_base = surface_idx + num_surf * 2 + 1;
    LeftXsurf_base_ext = surface_idx + num_surf * 3 + 1, RightXsurf_base_ext = surface_idx + num_surf * 4 + 1;
    LeftZsurf_base_ext = surface_idx + num_surf * 5 + 1, RightZsurf_base_ext = surface_idx + num_surf * 6 + 1;

    for (int i = 0; i < num_surf - 1; i++) {
        cmsh::model::geo::addLine(LeftXpoint[up][i + 1], RightXpoint[up][i + 1], line_idx + 1);

        cmsh::model::geo::addLine(LeftXpoint[up][i], LeftXpoint[up][i + 1], line_idx + 2);
        cmsh::model::geo::addLine(RightXpoint[up][i], RightXpoint[up][i + 1], line_idx + 3);

        cmsh::model::geo::addLine(LeftXpoint[down][i + 1], LeftXpoint[up][i + 1], line_idx + 4);
        cmsh::model::geo::addLine(RightXpoint[down][i + 1], RightXpoint[up][i + 1], line_idx + 5);

        Interface_line[NZ + 2].push_back(line_idx + 1);
        Interface_line[NZ + 2].push_back(line_idx + 2);
        Interface_line[NZ + 2].push_back(line_idx + 3);

        int Curve_tag5 = cmsh::model::geo::addCurveLoop({line_idx + 2, line_idx + 1, -(line_idx + 3), -LastyLine});
        int Curve_tag6 = cmsh::model::geo::addCurveLoop({LastzLine_1, line_idx + 2, -(line_idx + 4), -(LastxLine_1[i])});
        int Curve_tag7 = cmsh::model::geo::addCurveLoop({LastzLine_2, line_idx + 3, -(line_idx + 5), -(LastxLine_2[i])});

        LastxLine_1[i] = line_idx + 2;
        LastxLine_2[i] = line_idx + 3;

        // left_ext
        cmsh::model::geo::addLine(LeftXpoint[up][i + 1], LeftXpoint_ext[up][i + 2], line_idx + 6);

        cmsh::model::geo::addLine(LeftXpoint_ext[up][i + 1], LeftXpoint_ext[up][i + 2], line_idx + 7);
        cmsh::model::geo::addLine(LeftXpoint_ext[down][i + 2], LeftXpoint_ext[up][i + 2], line_idx + 8);
        // right_ext
        cmsh::model::geo::addLine(RightXpoint[up][i + 1], RightXpoint_ext[up][i + 2], line_idx + 9);

        cmsh::model::geo::addLine(RightXpoint_ext[up][i + 1], RightXpoint_ext[up][i + 2], line_idx + 10);
        cmsh::model::geo::addLine(RightXpoint_ext[down][i + 2], RightXpoint_ext[up][i + 2], line_idx + 11);

        Interface_line[NZ + 2].push_back(line_idx + 6);
        Interface_line[NZ + 2].push_back(line_idx + 7);

        Interface_line[NZ + 2].push_back(line_idx + 9);
        Interface_line[NZ + 2].push_back(line_idx + 10);

        Curve_tagl = cmsh::model::geo::addCurveLoop({LastLefty, line_idx + 7, -(line_idx + 6), -(line_idx + 2)});
        Curve_tagr = cmsh::model::geo::addCurveLoop({LastRighty, line_idx + 10, -(line_idx + 9), -(line_idx + 3)});

        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftZsurf_base_ext + i);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightZsurf_base_ext + i);

        Curve_tagl = cmsh::model::geo::addCurveLoop({LastzLine_3, line_idx + 7, -(line_idx + 8), -(LastxLine_1_ext[i + 1])});
        Curve_tagr = cmsh::model::geo::addCurveLoop({LastzLine_4, line_idx + 10, -(line_idx + 11), -(LastxLine_2_ext[i + 1])});

        LastxLine_1_ext[i + 1] = line_idx + 7;
        LastxLine_2_ext[i + 1] = line_idx + 10;

        cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftXsurf_base_ext + i);
        cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightXsurf_base_ext + i);

        cmsh::model::geo::addSurfaceFilling({Curve_tag5}, Zsurf_base + i);
        cmsh::model::geo::addSurfaceFilling({Curve_tag6}, LeftXsurf_base + i);
        cmsh::model::geo::addSurfaceFilling({Curve_tag7}, RightXsurf_base + i);

        surfaceTags.push_back(Zsurf_base + i);
        surfaceTags.push_back(LeftXsurf_base + i);
        surfaceTags.push_back(RightXsurf_base + i);

        LastSurface.push_back(Zsurf_base + i);
        LastLeftSurf.push_back(LeftZsurf_base_ext + i);
        LastRightSurf.push_back(RightZsurf_base_ext + i);

        LeftsurfaceTags.push_back(LeftXsurf_base + i);
        LeftsurfaceTags.push_back(LeftXsurf_base_ext + i);
        LeftsurfaceTags.push_back(LeftZsurf_base_ext + i);
        RightsurfaceTags.push_back(RightXsurf_base + i);
        RightsurfaceTags.push_back(RightXsurf_base_ext + i);
        RightsurfaceTags.push_back(RightZsurf_base_ext + i);

        LastyLine = line_idx + 1, LastLefty = line_idx + 6, LastRighty = line_idx + 9;
        LastzLine_1 = line_idx + 4;
        LastzLine_2 = line_idx + 5;
        LastzLine_3 = line_idx + 8;
        LastzLine_4 = line_idx + 11;
        line_idx += 11;
        surface_idx += 7;
    }

    cmsh::model::geo::addLine(LeftXpoint[up][num_surf - 1], LeftXpoint[up][num_surf], line_idx + 1);
    cmsh::model::geo::addLine(RightXpoint[up][num_surf - 1], RightXpoint[up][num_surf], line_idx + 2);

    Interface_line[NZ + 2].push_back(line_idx + 1);
    Interface_line[NZ + 2].push_back(line_idx + 2);

    Curve_tag8 = cmsh::model::geo::addCurveLoop({line_idx + 1, Line_idx6, -(line_idx + 2), -LastyLine});
    Curve_tag9 = cmsh::model::geo::addCurveLoop({LastzLine_1, line_idx + 1, -(Line_idx2), -(LastxLine_1[num_surf - 1])});
    Curve_tag10 = cmsh::model::geo::addCurveLoop({LastzLine_2, line_idx + 2, -(Line_idx4), -(LastxLine_2[num_surf - 1])});

    LastxLine_1[num_surf - 1] = line_idx + 1;
    LastxLine_2[num_surf - 1] = line_idx + 2;

    cmsh::model::geo::addSurfaceFilling({Curve_tag8} , Zsurf_base + num_surf - 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tag9} , LeftXsurf_base + num_surf - 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tag10}, RightXsurf_base + num_surf - 1);

    surfaceTags.push_back(Zsurf_base + num_surf - 1);
    surfaceTags.push_back(LeftXsurf_base + num_surf - 1);
    surfaceTags.push_back(RightXsurf_base + num_surf - 1);

    LastSurface.push_back(Zsurf_base + num_surf - 1);

    LeftsurfaceTags.push_back(LeftXsurf_base + num_surf - 1);
    RightsurfaceTags.push_back(RightXsurf_base + num_surf - 1);

    cmsh::model::geo::addLine(EndPoint_idx2, EndPoint_idx10, line_idx + 3);
    cmsh::model::geo::addLine(EndPoint_idx6, EndPoint_idx16, line_idx + 4);

    cmsh::model::geo::addLine(EndPoint_idx4, EndPoint_idx13, line_idx + 5);
    cmsh::model::geo::addLine(EndPoint_idx8, EndPoint_idx19, line_idx + 6);

    // left_ext
    cmsh::model::geo::addLine(LeftXpoint_ext[up][num_surf], LeftXpoint_ext[up][num_surf + 1], line_idx + 7);

    // right_ext
    cmsh::model::geo::addLine(RightXpoint_ext[up][num_surf], RightXpoint_ext[up][num_surf + 1], line_idx + 8);

    Interface_line[NZ + 2].push_back(line_idx + 3);
    Interface_line[NZ + 2].push_back(line_idx + 4);
    Interface_line[NZ + 2].push_back(line_idx + 5);
    Interface_line[NZ + 2].push_back(line_idx + 6);
    Interface_line[NZ + 2].push_back(line_idx + 7);
    Interface_line[NZ + 2].push_back(line_idx + 8);

    Curve_tagl = cmsh::model::geo::addCurveLoop({LastLefty, line_idx + 7, -Line_idx29, -(line_idx + 1)});
    Curve_tagr = cmsh::model::geo::addCurveLoop({LastRighty, line_idx + 8, -Line_idx32, -(line_idx + 2)});

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftZsurf_base_ext + num_surf - 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightZsurf_base_ext + num_surf - 1);

    Curve_tagl = cmsh::model::geo::addCurveLoop({LastzLine_3, line_idx + 7, -Line_idx10, -(LastxLine_1_ext[num_surf])});
    Curve_tagr = cmsh::model::geo::addCurveLoop({LastzLine_4, line_idx + 8, -Line_idx16, -(LastxLine_2_ext[num_surf])});

    LastxLine_1_ext[num_surf] = line_idx + 7;
    LastxLine_2_ext[num_surf] = line_idx + 8;

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, LeftXsurf_base_ext + num_surf - 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, RightXsurf_base_ext + num_surf - 1);

    LastLeftSurf.push_back(LeftZsurf_base_ext + num_surf - 1);
    LastRightSurf.push_back(RightZsurf_base_ext + num_surf - 1);

    LeftsurfaceTags.push_back(LeftXsurf_base_ext + num_surf - 1);
    LeftsurfaceTags.push_back(LeftZsurf_base_ext + num_surf - 1);

    RightsurfaceTags.push_back(RightXsurf_base_ext + num_surf - 1);
    RightsurfaceTags.push_back(RightZsurf_base_ext + num_surf - 1);

    // other surfaces
    surface_idx += 7;

    Curve_tag8 = cmsh::model::geo::addCurveLoop({line_idx + 3, Line_idx19, -(line_idx + 4), -Line_idx5});
    Curve_tag9 = cmsh::model::geo::addCurveLoop({Line_idx1, line_idx + 3, -Line_idx8, -LastLine_x1});
    Curve_tag10 = cmsh::model::geo::addCurveLoop({Line_idx3, line_idx + 4, -Line_idx14, -LastLine_x2});

    cmsh::model::geo::addSurfaceFilling({Curve_tag8}, surface_idx + 1);
    cmsh::model::geo::addSurfaceFilling({Curve_tag9}, surface_idx + 2);
    cmsh::model::geo::addSurfaceFilling({Curve_tag10}, surface_idx + 3);

    FrontsurfaceTags.push_back(surface_idx + 1);
    FrontsurfaceTags.push_back(surface_idx + 2);
    FrontsurfaceTags.push_back(surface_idx + 3);

    LeftsurfaceTags.push_back(surface_idx + 2);
    RightsurfaceTags.push_back(surface_idx + 3);

    Curve_tag8 = cmsh::model::geo::addCurveLoop({line_idx + 5, Line_idx20, -(line_idx + 6), -Line_idx6});
    Curve_tag9 = cmsh::model::geo::addCurveLoop({Line_idx2, line_idx + 5, -Line_idx11, -LastLine_x3});
    Curve_tag10 = cmsh::model::geo::addCurveLoop({Line_idx4, line_idx + 6, -Line_idx17, -LastLine_x4});

    cmsh::model::geo::addSurfaceFilling({Curve_tag8}, surface_idx + 4);
    cmsh::model::geo::addSurfaceFilling({Curve_tag9}, surface_idx + 5);
    cmsh::model::geo::addSurfaceFilling({Curve_tag10}, surface_idx + 6);

    BacksurfaceTags.push_back(surface_idx + 4);
    BacksurfaceTags.push_back(surface_idx + 5);
    BacksurfaceTags.push_back(surface_idx + 6);

    LeftsurfaceTags.push_back(surface_idx + 5);
    RightsurfaceTags.push_back(surface_idx + 6);

    Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx29, Line_idx27, -Line_idx28, -(line_idx + 5)});
    Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx32, Line_idx30, -Line_idx31, -(line_idx + 6)});

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 7);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 8);

    Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx23, Line_idx21, -Line_idx22, -(line_idx + 3)});
    Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx26, Line_idx24, -Line_idx25, -(line_idx + 4)});

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 9);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 10);

    LeftsurfaceTags.push_back(surface_idx + 7);
    LeftsurfaceTags.push_back(surface_idx + 9);
    LastLeftSurf.push_back(surface_idx + 7);
    LastLeftSurf.push_back(surface_idx + 9);

    RightsurfaceTags.push_back(surface_idx + 8);
    RightsurfaceTags.push_back(surface_idx + 10);
    LastRightSurf.push_back(surface_idx + 8);
    LastRightSurf.push_back(surface_idx + 10);

    Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx10, Line_idx27, -Line_idx12, -(LastxLine_1_ext[num_surf + 1])});
    Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx16, Line_idx30, -Line_idx18, -(LastxLine_2_ext[num_surf + 1])});

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 11);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 12);

    Curve_tagl = cmsh::model::geo::addCurveLoop({Line_idx7, Line_idx21, -Line_idx9, -(LastxLine_1_ext[0])});
    Curve_tagr = cmsh::model::geo::addCurveLoop({Line_idx13, Line_idx24, -Line_idx15, -(LastxLine_2_ext[0])});

    cmsh::model::geo::addSurfaceFilling({Curve_tagl}, surface_idx + 13);
    cmsh::model::geo::addSurfaceFilling({Curve_tagr}, surface_idx + 14);

    LeftsurfaceTags.push_back(surface_idx + 11);
    LeftsurfaceTags.push_back(surface_idx + 13);
    RightsurfaceTags.push_back(surface_idx + 12);
    RightsurfaceTags.push_back(surface_idx + 14);

    LastLine_x1 = line_idx + 3, LastLine_x2 = line_idx + 4;
    LastLine_x3 = line_idx + 5, LastLine_x4 = line_idx + 6;
    LastFrontSurf = surface_idx + 1;
    LastBackSurf = surface_idx + 4;

    LastxLine_1_ext[0] = Line_idx21;
    LastxLine_1_ext[num_surf + 1] = Line_idx27;
    LastxLine_2_ext[0] = Line_idx24;
    LastxLine_2_ext[num_surf + 1] = Line_idx30;

    line_idx += 8;
    surface_idx += 14;

    cmsh::model::geo::addSurfaceLoop(surfaceTags, volume_idx + 1);
    cmsh::model::geo::addVolume({volume_idx + 1}, volume_idx + 1);

    cmsh::model::geo::addSurfaceLoop(FrontsurfaceTags, volume_idx + 2);
    cmsh::model::geo::addVolume({volume_idx + 2}, volume_idx + 2);

    cmsh::model::geo::addSurfaceLoop(BacksurfaceTags, volume_idx + 3);
    cmsh::model::geo::addVolume({volume_idx + 3}, volume_idx + 3);

    cmsh::model::geo::addSurfaceLoop(LeftsurfaceTags, volume_idx + 4);
    cmsh::model::geo::addVolume({volume_idx + 4}, volume_idx + 4);

    cmsh::model::geo::addSurfaceLoop(RightsurfaceTags, volume_idx + 5);
    cmsh::model::geo::addVolume({volume_idx + 5}, volume_idx + 5);

    cellzone_id[2].push_back(volume_idx + 1);
    cellzone_id[3].push_back(volume_idx + 2);
    cellzone_id[3].push_back(volume_idx + 3);
    cellzone_id[3].push_back(volume_idx + 4);
    cellzone_id[3].push_back(volume_idx + 5);

    volume_idx += 5;
    LastPoint_1 = EndPoint_idx2;
    LastPoint_2 = EndPoint_idx4;
    LastPoint_3 = EndPoint_idx6;
    LastPoint_4 = EndPoint_idx8;
    LastLine_1 = Line_idx5;
    LastLine_2 = Line_idx6;

    EnterLastPoint.push_back(point_idx);

    Interface_surf[NZ + 2].insert(Interface_surf[NZ + 2].end(), LastSurface.begin(), LastSurface.end());
    Interface_surf[NZ + 2].insert(Interface_surf[NZ + 2].end(), LastLeftSurf.begin(), LastLeftSurf.end());
    Interface_surf[NZ + 2].insert(Interface_surf[NZ + 2].end(), LastRightSurf.begin(), LastRightSurf.end());
    Interface_surf[NZ + 2].push_back(LastFrontSurf);
    Interface_surf[NZ + 2].push_back(LastBackSurf);

    cmsh::model::geo::synchronize();

    cmsh::model::mesh::generate(3);

    /* 4、Output */
    std::vector<GridPoint<datatype>> dstpointArray;
    std::vector<GridCell<int, datatype, 4>> dstcellArray;

    std::vector<std::size_t> nodeTags;
    std::vector<double> nodeCoords, parametricCoords;
    cmsh::model::mesh::getNodes(nodeTags, nodeCoords, parametricCoords);
    
    dstpointArray.resize(nodeTags.size());
#pragma omp parallel
    for(std::size_t i = 0; i < nodeTags.size(); ++i) {
        dstpointArray[i] = GridPoint<datatype>(nodeCoords[i * 3], nodeCoords[i * 3 + 1], nodeCoords[i * 3 + 2]);
    }

    /* Origin Point */
    cmsh::model::mesh::getNodes(nodeTags, nodeCoords, parametricCoords, 0);
    int idx = 1;

    for (int i = 0 ; i < EnterLastPoint.size(); i++) {
       for (; idx <= EnterLastPoint[i]; idx++) {
            dstpointArray[idx - 1]._interfaceIdx = 11 + NZ - i;
       }
    }

    /* Line Point */
    for (int i = 0; i < NZ + 3; i ++)
    {   
        for (int j = 0 ; j < Interface_line[i].size(); j++)
        {
            nodeTags.clear();
            cmsh::model::mesh::getNodes(nodeTags, nodeCoords, parametricCoords, 1, Interface_line[i][j]);
            // std::cout << "Line " << Interface_line[i][j] << " has " << nodeTags.size() << " nodes." << std::endl;
            for (int k = 0; k < nodeTags.size(); k++) {
                idx = nodeTags[k] - 1;
                dstpointArray[idx]._interfaceIdx = 11 + NZ - i;
            }
        }
    }

    /* Surface Point */
    for (int i = 0; i < NZ + 3; i ++)
    {
        for (int j = 0 ; j < Interface_surf[i].size(); j++)
        {
            nodeTags.clear();
            // printf("tmp_surf_idx: %d\n", tmp_surf_idx);
            cmsh::model::mesh::getNodes(nodeTags, nodeCoords, parametricCoords, 2, Interface_surf[i][j]);
            // std::cout << "Surface " << Interface_surf[i][j] << " has " << nodeTags.size() << " nodes." << std::endl;
            for (int k = 0; k < nodeTags.size(); k++) {
                idx = nodeTags[k] - 1;
                dstpointArray[idx]._interfaceIdx = 11 + NZ - i;
            }
        }
    }

    /* Fault Point */
    for (int i = 0; i < FaultArray.size(); i++) {
        nodeTags.clear();
        cmsh::model::mesh::getNodes(nodeTags, nodeCoords, parametricCoords, 2, FaultSurface_idx[i]);
        std::cout << "Fault " << i << " has " << nodeTags.size() << " nodes." << std::endl;
        for (int j = 0; j < nodeTags.size(); j++) {
            idx = nodeTags[j] - 1;
            dstpointArray[idx]._interfaceIdx = i + 1;
        }
    }

    dstgrid.SetPointArray(dstpointArray);

    std::vector<std::size_t> elementTags;
    std::vector<std::size_t> cellnodeTags;
    cmsh::model::mesh::getElementsByType(4, elementTags, cellnodeTags);

    dstcellArray.resize(elementTags.size());

    int tmp_offset = 0;
    for(int i = 0; i < 6; i++) {
        for(int j = 0; j < cellzone_id[i].size(); j++) {
                elementTags.clear();
                cellnodeTags.clear();
                cmsh::model::mesh::getElementsByType(4, elementTags, cellnodeTags, cellzone_id[i][j]);

#pragma omp parallel
                for (int k = 0; k < elementTags.size(); k++) {
                    std::array<int, 4> tmp_points = {
                        static_cast<int>(cellnodeTags[k * 4]) - 1, 
                        static_cast<int>(cellnodeTags[k * 4 + 1]) - 1, 
                        static_cast<int>(cellnodeTags[k * 4 + 2]) - 1, 
                        static_cast<int>(cellnodeTags[k * 4 + 3]) - 1
                    };
                    dstcellArray[tmp_offset + k] = GridCell<int, datatype, 4>((tmp_points), 0.0, 0.0, 0.0, 0.0, i + 1);
                }
                tmp_offset += elementTags.size();
        }

        if(i == 0)
        {
            dstgrid.SetNumOriginCell(tmp_offset);
        }
    }

    dstgrid.SetCellArray(dstcellArray);

    cmsh::finalize();

    return 0;   
}

template int GridRemesh<double>(const CornerGrid<false, double, GridCell<GridPoint<double>, double, 8>>&, 
                                FinElemGrid<GridCell<int, double, 4>, double, 4>&, 
                                double, double, double, double, double);

template int GridRemesh<float>(const CornerGrid<false, float, GridCell<GridPoint<float>, float, 8>>&, 
                                FinElemGrid<GridCell<int, float, 4>, float, 4>&, 
                                double, double, double, double, double);

template <typename datatype>
int ComputePrefix(
    const CornerGrid<true, datatype, GridCell<GridPoint<datatype>, datatype, 8>>& finegrid, 
    int NX, 
    int NY, 
    int NZ, 
    std::vector<std::pair<datatype, datatype>>& prefixZ, 
    std::vector<std::vector<std::pair<datatype, datatype>>>& prefixY) 
{
    auto src_CellArray = finegrid.GetCellArray();

    prefixZ.resize(NZ);
    prefixY.resize(NZ);

    // Assume x-y-z Cycle
    for (int k = 0; k < NZ; k++) {
        datatype z_max = -1e10, z_min = 1e10;
        prefixY[k].resize(NY);
        for (int j = 0; j < NY; j++) {
            datatype y_max = -1e10, y_min = 1e10;
            for (int i = 0; i < NX; i++) {
                int cellIndex = k * NY * NX + j * NX + i; // Compute cell index
                for (int index = 0; index < 8; index++) // Process 8 corners of the cell
                {
                    datatype py = src_CellArray[cellIndex]._points[index]._y;
                    datatype pz = src_CellArray[cellIndex]._points[index]._z;
                    if (pz > z_max) z_max = pz;
                    if (pz < z_min) z_min = pz;
                    if (py > y_max) y_max = py;
                    if (py < y_min) y_min = py;
                }
            }
            prefixY[k][j] = std::make_pair(y_min, y_max);
        }
        prefixZ[k] = std::make_pair(z_min, z_max);
    }

    return 0;
}

template int ComputePrefix<double>(const CornerGrid<true, double, GridCell<GridPoint<double>, double, 8>>&, 
                                    int, int, int, 
                                    std::vector<std::pair<double, double>>&, 
                                    std::vector<std::vector<std::pair<double, double>>>&);

template int ComputePrefix<float>(const CornerGrid<true, float, GridCell<GridPoint<float>, float, 8>>&, 
                                    int, int, int, 
                                    std::vector<std::pair<float, float>>&, 
                                    std::vector<std::vector<std::pair<float, float>>>&);

template <typename datatype>
static inline void _SearchMap(GridCell<int, datatype, 4>& cell, 
                              int NX, int NY, int NZ,
                              datatype center_x, datatype center_y, datatype center_z,
                              const std::vector<GridCell<GridPoint<datatype>, datatype, 8>>& src_CellArray,
                              const std::vector<std::pair<datatype, datatype>>& prefixZ, 
                              const std::vector<std::vector<std::pair<datatype, datatype>>>& prefixY) 
{   
    std::vector<int> z_idx;
    for (int i = 0; i < prefixZ.size(); i++) {
        if (center_z >= prefixZ[i].first && center_z <= prefixZ[i].second) {
            z_idx.push_back(i);
        }
    }

    for (int i = 0; i < z_idx.size(); i++) {
        int z = z_idx[i];
        // std::cout << "z: " << z << std::endl;
        std::vector<int> y_idx;
        for (int j = 0; j < prefixY[z].size(); j++) {
            if (center_y >= prefixY[z][j].first && center_y <= prefixY[z][j].second) {
                y_idx.push_back(j);
            }
        }

        for (int j = 0; j < y_idx.size(); j++) {
            int y = y_idx[j];
            // std::cout << "y: " << y << std::endl;
            for (int k = 0; k < NX; k++) {
                int cellIndex = z * NY * NX + y * NX + k;
                auto src_cell = src_CellArray[cellIndex];
                datatype x_max = -1e10, x_min = 1e10;
                datatype y_max = -1e10, y_min = 1e10;
                datatype z_max = -1e10, z_min = 1e10;

                for (int index = 0; index < 8; index++) {
                    datatype px = src_cell._points[index]._x;
                    datatype py = src_cell._points[index]._y;
                    datatype pz = src_cell._points[index]._z;
                    if (px > x_max) x_max = px;
                    if (px < x_min) x_min = px;
                    if (py > y_max) y_max = py;
                    if (py < y_min) y_min = py;
                    if (pz > z_max) z_max = pz;
                    if (pz < z_min) z_min = pz;
                }
                if (center_x >= x_min && center_x <= x_max
                    && center_y >= y_min && center_y <= y_max
                    && center_z >= z_min && center_z <= z_max) {
                    cell._youngMod = src_cell._youngMod;
                    cell._poissonRatio = src_cell._poissonRatio;
                    cell._rockDensity = src_cell._rockDensity;
                    cell._porosity = src_cell._porosity;

                    // std::cout << "cell._youngMod: " << cell._youngMod << " ,src_cell._youngMod: " << src_cell._youngMod << std::endl;
                    // std::cout << "cell._poissonRatio: " << cell._poissonRatio << " ,src_cell._poissonRatio: " << src_cell._poissonRatio << std::endl;
                    // std::cout << "cell._rockDensity: " << cell._rockDensity << " ,src_cell._rockDensity: " << src_cell._rockDensity << std::endl;
                    // std::cout << "cell._porosity: " << cell._porosity << " ,src_cell._porosity: " << src_cell._porosity << std::endl;
                    return;
                }
            }
        }
    }
    
}

template <typename datatype>
int PropMapping(const CornerGrid<true, datatype, GridCell<GridPoint<datatype>, datatype, 8>>& srcgrid,
                const std::vector<std::pair<datatype, datatype>>& prefixZ, 
                const std::vector<std::vector<std::pair<datatype, datatype>>>& prefixY,
                FinElemGrid<GridCell<int, datatype, 4>, datatype, 4>& dstgrid)
{
    const auto& src_CellArray = srcgrid.GetCellArray();

    auto dst_CellArray = dstgrid.GetCellArray();
    const auto& dst_PointArray = dstgrid.GetPointArray();

    /* Origin Grid Mapping */
#pragma omp parallel
    for (int i = 0; i < dstgrid.GetNumOriginCell(); i++) {
        datatype center_x = 0.0, center_y = 0.0, center_z = 0.0;
        for (int j = 0; j < 4; j++) {
            int idx = dst_CellArray[i]._points[j];
            datatype px = dst_PointArray[idx]._x;
            datatype py = dst_PointArray[idx]._y;
            datatype pz = dst_PointArray[idx]._z;

            center_x += px;
            center_y += py;
            center_z += pz;
        }
        center_x /= 4;
        center_y /= 4;
        center_z /= 4;

        _SearchMap<datatype>(dst_CellArray[i], srcgrid.GetNx(), srcgrid.GetNy(), srcgrid.GetNz(), center_x, center_y, center_z, src_CellArray, prefixZ, prefixY);
    }

    /* Extension Grid Mapping */
    // std::cout << "Total Cell is " << dst_CellArray.size() << std::endl;
    // std::cout << "Oversides offset is " << dstgrid.GetNumOriginCell() << std::endl;
    // std::cout << "Overburden offset is " << dstgrid.GetOverburdenOffset() << std::endl;
    // std::cout << "Undersides offset is " << dstgrid.GetUndersidesOffset() << std::endl;
    // std::cout << "Underburden offset is " << dstgrid.GetUnderburdenOffset() << std::endl;
    // std::cout << "Sidesburden offset is " << dstgrid.GetSidesburdenOffset() << std::endl;

    std::vector<KDNode<datatype>> gridDataArray;
    gridDataArray.resize(src_CellArray.size());

#pragma omp parallel for
    for (int i = 0; i < src_CellArray.size(); i++) {
        datatype center_x = 0.0, center_y = 0.0, center_z = 0.0;
        for (int j = 0; j < 8; j++) {
            datatype px = src_CellArray[i]._points[j]._x;
            datatype py = src_CellArray[i]._points[j]._y;
            datatype pz = src_CellArray[i]._points[j]._z;

            center_x += px;
            center_y += py;
            center_z += pz;
        }
        center_x /= 8;
        center_y /= 8;
        center_z /= 8;

        gridDataArray[i] = KDNode<datatype>(center_x, center_y, center_z, src_CellArray[i]._youngMod, src_CellArray[i]._poissonRatio, src_CellArray[i]._rockDensity, src_CellArray[i]._porosity);
    }

    KDTree3D<KDNode<datatype>, datatype> tree(gridDataArray);

#pragma omp parallel for
    for (int i = dstgrid.GetNumOriginCell(); i < dst_CellArray.size(); i++) {
        datatype center_x = 0.0, center_y = 0.0, center_z = 0.0;
        for (int j = 0; j < 4; j++) {
            int idx = dst_CellArray[i]._points[j];
            datatype px = dst_PointArray[idx]._x;
            datatype py = dst_PointArray[idx]._y;
            datatype pz = dst_PointArray[idx]._z;

            center_x += px;
            center_y += py;
            center_z += pz;
        }
        center_x /= 4;
        center_y /= 4;
        center_z /= 4;

        KDNode<datatype> tmp(center_x, center_y, center_z, 0.0, 0.0, 0.0, 0.0);
        tree.K_Interpolation(tmp, 8);
        dst_CellArray[i]._youngMod = tmp._youngMod;
        dst_CellArray[i]._poissonRatio = tmp._poissonRatio;
        dst_CellArray[i]._rockDensity = tmp._rockDensity;
        dst_CellArray[i]._porosity = tmp._porosity;
    // std::cout << "i: " << i << " ,youngMod: " << dst_CellArray[i]._youngMod << " ,poissonRatio: " << dst_CellArray[i]._poissonRatio << " ,rockDensity: " << dst_CellArray[i]._rockDensity << " ,porosity: " << dst_CellArray[i]._porosity << std::endl;
    }

    dstgrid.SetCellArray(dst_CellArray);

    return 0;
}

template int PropMapping<double>(const CornerGrid<true, double, GridCell<GridPoint<double>, double, 8>>&,
                                 const std::vector<std::pair<double, double>>&,
                                 const std::vector<std::vector<std::pair<double, double>>>&,
                                 FinElemGrid<GridCell<int, double, 4>, double, 4>&);

template int PropMapping<float>(const CornerGrid<true, float, GridCell<GridPoint<float>, float, 8>>&, 
                                const std::vector<std::pair<float, float>>&,
                                const std::vector<std::vector<std::pair<float, float>>>&,
                                FinElemGrid<GridCell<int, float, 4>, float, 4>&);